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Computational  Support  for  the  Study  of  Lifetime  Distribution 

Characteristics 

R.  R.  Read 
ABSTRACT 

The  report  develops  and  makes  available  programs 
that  treat  the  support  functions  of  a  set  of  survivor 
distributions:  Weibull,  Gamma,  and  Lognormal.  The  issues 
of  model  characterization  functions,  maximum  likelihood 
estimation,  bias  reduction,  and  censored  samples  are  treated 
generally.  The  general  material  is  made  explicit  for  the 
distributions  named.  It  features  open  code,  allowing  the 
user  to  pursue  plans  of  his  own. 

The  paper  also  contains  some  items  of  broader 
interest.  First,  a  technique  is  developed  that  offers 
substantial  reduction  in  the  dependence  of  the  initialization 
values  for  the  success  of  the  Newton-Raphson  iteration 
technique.  Second,  high-precision,  numerical  analysis 
techniques  are  developed  for  the  parallel  computation  of 
several  derivatives  of  the  Incomplete  Gamma  function. 


1.  Introduction 


Much  has  been  written  on  life  length  distribution  applieations — espeeially  those 
in  the  field  of  reliability.  The  first  eourse  textbooks  usually  eover  the  ground  well,  but  the 
examples  often  are  shallow  beeause  there  are  so  few  distributions  that  ean  readily  be 
managed,  supported  only  by  the  ealeulus.  The  eomputational  aspeets  are  extremely  easy 
when  dealing  with  the  exponential  distribution,  and  they  beeome  quite  diffieult  when 
using  other  eommon  distributions.  This  situation  ereates  difficulties  in  the  teaching  of 
first  courses,  the  writing  of  textbooks,  and,  of  eourse,  in  the  execution  of  research.  The 
use  of  the  exponential  distribution  beeomes  overworked  and  dull,  and  its  sole  use  masks 
the  dangers  that  ean  be  eneountered. 

The  present  report  offers  useful  relief  to  these  problems.  Algorithms  and  S-Plus 
code  are  developed  so  that  the  user  ean  explore  on  the  eomputer  the  eonsequenees  of 
using  the  Weibull,  Gamma,  and  Lognormal  distributions.  Also,  a  general  mathematieal 
structure  is  presented  so  that  a  programmer  ean  extend  and  deal  with  other  distribution 
families. 

At  the  same  time,  we  are  beginning  to  see  produets  of  this  type  from  other 
software  writers.  The  Splida  system  by  Meeker  and  Alvarez  has  eompleted  beta  testing 
and  has  appeared  on  the  market.  This  system  appears  to  be  very  eomprehensive.  It  is 
written  in  S-Plus  and  features  a  large  number  of  speeialized  drop-down  dialog  menus.  A 
produet  named  Reliasoft  is  already  on  the  market.  It  too  requires  that  the  user  fill  out 
dialog  boxes.  The  present  work  is  far  less  ambitious,  but  eontains  open  S-Plus  eode  that 
allows  the  user  to  deal  direetly  with  issues  of  his  own  ehoosing  and  to  eheek  the  preeision 
of  his  solutions.  It  eontains  a  number  of  eommand  line  funetions  that  a  user  ean  integrate 
into  his  own  specialized  problem-solving  paekage.  (Meeker  offers  to  supply  computer 
eode  upon  request;  it  is  not  known  whether  this  ineludes  the  supporting  algorithmic 
analysis.) 

Attention  is  restricted  to  the  eontinuous  time  life  length  variables.  General 
formulae  are  developed  for  these.  A  first  goal  is  to  eompute  and  graph  five  basic  support 
functions,  namely. 


f(t)  =  probability  density  functions  (pdf), 

S(t)  =  survivor  function, 
h(t)  =  hazard  function, 

H(t)  =  cumulative  hazard  funetion, 

LL(t)=  expeeted  residual  life. 

The  first  two  often  are  ineluded  in  the  standard  statistieal  software  paekages.  The 
others  usually  require  the  praetitioner  to  exereise  some  ealeulus  and  programming.  The 
user  should  be  familiar  with  the  shapes  of  these  functions  and  how  they  vary  from  model 
to  model.  In  fact,  model  seleetion  is  often  a  judgment  call  rather  than  a  statistieal 
exereise.  In  terms  of  these  support  funetions,  some  model  eomparative  graphing  ean  be 
quite  valuable. 
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Complete  samples  are  those  for  whieh  the  testing  proeess  is  not  finished  until  all 
objeets  have  realized  their  entire  lifetimes.  In  these  eases,  the  stoek  statistieal  methods  are 
utilized.  One  usually  ealls  upon  maximum  likelihood  to  serve  for  parameter  estimation. 
These  have  well  studied  properties.  They  are  not  easy  to  find  and  often  iterative  methods 
are  created  for  this  purpose.  Seldom  are  these  available  in  the  standard  statistical  software 
packages.  The  report  contains  examples  of  how  one  can  generate  such  methods.  The 
large  sample  theory  for  maximum  likelihood  estimators  is  well  established  and  it  is  used 
to  create  confidence  regions  for  the  parameters. 

Generally,  maximum  likelihood  estimators  are  known  to  be  biased.  This  issue  has 
not  received  great  attention.  The  methods  for  dealing  with  it  are  difficult  and  the  effects 
are  not  believed  to  be  great,  except  possibly  for  small  sample  sizes.  But  the  report  does 
consider  a  first  order  bias  reduction  method  and  applies  it  to  the  models  treated.  At  least 
the  reader  is  allowed  to  view  the  nature  of  the  problem. 

Recent  times  have  seen  the  development  of  methods  to  manage  censored  life 
length  data,  i.e.,  data  collection  plans  that  do  not  collect  full  life  length  information  on  all 
subjects.  We  consider  right-censored  data  only,  i.e.,  each  subject  is  either  observed  to 
expire  or  has  survived  beyond  a  known  point.  Plans  that  await  a  complete  data  set  are 
often  too  expensive.  The  maximum  likelihood  point  estimation  schemes  are  available  in 
concept,  but  the  implementation  is  often  intricate  and  requires  special  methods.  The 
iteration  functions  utilized  are  advertised  as  requiring  high  quality  initialization  points.  A 
method  is  proposed  that  appears  to  give  substantial  relief  to  this  problem  and  shows  great 
promise  at  this  point  in  time. 

When  dealing  with  the  gamma  distribution  in  the  censored  case,  one  must  face  the 
computation  of  the  derivatives  of  the  incomplete  gamma  function.  Methods  for 
accomplishing  this  are  developed. 

The  report  is  organized  as  follows.  Section  2  contains  the  general  mathematical 
structure  for  the  support  functions,  for  complete  samples,  for  a  bias  reduction  method  and 
for  censored  samples.  This  is  followed,  in  Section  3,  by  the  presentation  of  some 
computations  and  graphs  typical  of  the  capabilities  that  are  being  supported.  This  part 
does  much  to  illustrate  the  use  of  the  programs  and  the  kinds  of  issues  that  may  be 
addressed.  The  reader  who  is  well  versed  in  the  issues  will  find  this  chapter  most  useful 
as  it  identifies  the  programs  and  illustrates  their  use.  Sections  4,  5,  and  6  treat  the 
mathematical  details  for  the  three  popular  distributions:  Weibull,  Gamma,  and 
Lognormal,  respectively.  Explicit  formulae  are  developed  for  the  characterizing  support 
functions,  the  treatment  of  complete  samples  with  maximum  likelihood,  bias  reduction, 
and  methods  for  censored  samples. 

Two  data  sets  have  been  enlisted  to  test  the  programs  in  Section  3.  The  Lieblein 
and  Zelen  (1956)  ball-bearing  data  is  utilized  for  complete  samples.  It  consists  of  23 
failure  times  of  ball  bearings  measured  in  millions  of  revolutions.  It  has  been  exploited 
by  [Meeker  and  Escobar,  p.  4]  to  illustrate  the  model-fitting  problem.  A  similar  use  is 
presented  in  Section  3.  The  data  are 
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BallB;  17.88  28.92  33.00  41.52  42.12  45.60  48.40  51.84  51.96  54.12  67.80 
68.64  68.64  68.88  84.12  93.12  98.64  105.12  105.84  127.92  128.04  173.40 


A  second  data  set  was  chosen  to  test  the  methods  that  involve  censoring.  It  is  the 
set  of  survival  times  for  multiple  myeloma  patients  from  [Lawless,  p.  337]  consisting  of 
48  complete  lifetimes,  x  and  17  right-censored  values,  t.  Both  are  measured  in  months. 
The  data  are 

myeloma: 

x:  1  1  2  2  2  3  5  5  6  6  6  6  7  7  7  9  11  11  11  11  11  13  14 
15  16  16  17  17  18  19  19  24  25  26  32  35  37  41  42  51  52  54  58  66  67  88 
89  92 

t:  4  4  7  7  8  12  11  12  13  16  19  19  28  41  53  57  77 

Fuller  details  of  implementation  of  the  methods  are  contained  in  five  appendices. 

They  treat  specialized  subjects  and  each  contains  supporting  S-Plus  listings.  Their 
contents  are: 

Appendix  A.  Asymptotic  Expansions  for  the  Polygamma  Functions. 

Appendix  B.  Analysis  for  Computational  Support  of  the  Weibull  Distribution. 
Appendix  C.  Implementation  of  the  General  Censored  Data  Estimation  Scheme. 
Appendix  D.  Derivatives  of  the  Incomplete  Gamma  Function. 

Appendix  E.  S-Plus  Eistings  of  Miscellaneous  Code. 

The  programs  can  be  obtained  in  electronic  form  by  contacting  the  author. 
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2,  General  Formula 


The  specialized  programs  follow  a  general  structure,  which  is  presented  here.  A 
number  of  computational  situations  are  envisioned.  Only  the  continuous  case  is  treated. 
First,  we  deal  with  important  support  functions  used  in  reliability  and  lifetime  analysis. 
This  includes  a  visit  to  the  optimal  planned  replacement  policy  formulae.  Second,  the 
basic  maximum  likelihood  estimation  procedures  for  complete  samples  is  presented.  This 
includes  iteration  schemes  for  finding  the  solution,  the  information  matrix,  and  the 
construction  of  asymptotic  confidence  regions.  This  also  includes  the  method  for 
reducing  the  sensitivity  of  the  initialization  values.  Third,  the  development  of  a  general 
technique  for  bias  reduction  is  considered.  Fourth,  the  general  structure  for  dealing  with 
maximum  likelihood  estimation  for  right-censored  samples  is  presented. 

Code  is  needed  for  the  exploration  of  the  properties  of  various  models. 
Relationships  among  the  various  functions  are  summarized  [Leemis,  p.  55].  Those  that 
are  most  useful  in  the  present  work  are  listed.  A  continuous  life  length  random  variable  X 
has  a  probability  density  function,  f(t)  and  a  cumulative  distribution  function  F(t).  From 


these  one  can  characterize  others. 

2a.  Model  Characterization  Functions 

The  survivor  function: 

S(t)  =  f(u)du  =  l-F(t)  =  Pr{X>t}.  (2.1) 

The  hazard  function,  age  specific  failure  rate: 

h(t)  =  f(t)/S(t)  =  H'(t).  (2.2) 

The  cumulative  hazard  function: 

H(t)  =  [‘h(u)du  =  -ln[S(t)].  (2.3) 

J  0 

The  integrated  survivor  function: 

C  00 

SS(t)  =  S(u)du.  (2.4) 

The  mean  residual  life: 

C  00 

LL(t)  =  uf(u)du/S(t)  =  SS(t)/S(t)  =  E{X-t|X>t}.  (2.5) 

2b,  Cost  Calculation  of  Planned  Replacement  Policies 


Let  us  examine  the  effect  of  model  choice  when  a  planned  replacement  policy  is 
considered.  Such  policies  can  be  advantageous  when  dealing  with  an  IFR  system  and  the 
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cost  of  a  planned  replacement  is  lower  than  the  cost  of  replacing  at  failure.  The  structure 
of  the  calculation  is  developed  in  [Barlow  and  Proschan], 

Let  k  >  0  be  the  cost  of  a  planned  replacement,  but  if  the  replacement  must 
happen  because  of  a  device  failure  then  an  additional  cost  c  is  incurred  to  be  added  to  k. 
Let  Rn  be  the  cost  of  the  n*’’  replacement,  and  Xn  is  the  lifelength  of  the  n*  unit  placed  in 
service.  The  system  [Xn]  is  assumed  to  be  i.i.d.  The  cost  is  random  because  it  depends 
on  the  type  of  replacement. 

Rn  =  k  if  Xn  >  t  and  Rn  =  k  +  c  if  Xn  <  t.  (2.6) 

Now  we  can  describe  the  cost  for  the  interval  (0,  s]. 

This  cost  is  a  random  variable,  but  the  long-term  average  cost  stabilizes.  That  is 

Z{s)  E{R,} 

^  as  s  ^  00. 

5  E{X,} 

The  main  formula  for  long-term  average  cost  using  replacement  policy  t,  i.e.,  replace  at 
age  t,  is 


C(t)  = 


k  +  cF(t) 

f  ^  ^ 

^^S{w)dw 


where  S(w)  is  the  survivor  function  of  the  distribution 
expression  may  also  be  computed  using  p,  -  SS(t)  and  p  = 
and  Proschan,  1981]  or  [Prentice  and  Kalbfleisch]. 


(2.7) 

F.  The  denominator  of  this 
E{X}  =  SS(0).  See  [Barlow 


2c.  The  Likelihood  Equations  and  Solution  Technique;  Structural  Overview 

The  structure  of  the  Itkelihood  function  has  different  appearances  depending  upon 
whether  the  data  are  complete  or  censored.  There  are  a  few  general  points  that  apply  to 
either  case.  These  are  presented  first,  followed  by  representations  of  how  the  details  can 
change  depending  upon  the  two  cases.  The  parameter  0  may  be  a  vector  of  several 
components.  Let  f  (0)  be  the  likelihood  function,  and 


L(0)  =  log[f(0)]  (2.8) 

be  its  logarithm.  We  are  concerned  with  the  smooth  settings  in  which  the  maximum 
likelihood  estimates  are  found  using  gradient  methods.  The  partial  derivative  with  respect 
to  the  individual  members  of  0,  when  the  particular  subscript(s)  plays  no  immediate  role 
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will  be  marked  as  Le;  Lee;  and  L  ;  the  latter  case  representing  any  pair  of  the  mixed 
partial  derivatives. 

The  partial  derivative  of  L  with  respect  to  0  has  expected  value  equal  to  zero  and 
the  expected  square  of  this  quantity  is  the  negative  of  the  expected  value  of  the  second 
partial  of  the  log  likelihood.  The  requisite  smoothness  conditions  for  interchanging  the 
expectation  operation  with  the  appropriate  partial  derivatives  are  presumed. 

The  matrix  {L  gg  }  is  known  as  the  Hessian  and  the  negative  of  its  expected 

value  is  the  (Fisher)  information  matrix.  Call  the  former  H,  the  latter  nio  (use  lo  for  the 
information  in  a  single  observation),  and  Le  the  vector  of  partials  of  the  log  likelihood. 
Often  the  system  of  equations  Le  =  0  must  be  solved  by  iterative  methods.  This  is  done 
using  a  Newton-Raphson  type  technique.  Two  options  are  presented; 

0{ku)  ^  0(k)  _  pj-1  and  +  (nIo)-'  Lg .  (2.9) 

Termination  of  the  iteration  occurs  when  there  is  no  change  in  the  maximum  value.  Both 
methods  require  a  good  initialization,  0*^°\  To  some,  there  appears  to  be  empirical 
evidence  to  use  the  second  choice,  when  feasible,  but  the  calculation  of  lo  is  often 
difficult. 

Since  convergence  of  the  iteration  function  is  sensitive  to  the  initialization,  0*^^*^ 
the  author  found  relief  from  this  problem  by  utilizing  a  golden  section  search  along  the 
segment  (0  0  Typically,  the  iteration  steps  over-swing  the  maximum  in  this 

direction,  often  by  a  factor  of  two  and  sometimes  by  a  factor  of  10.  It  can  pay  to  seek  a 

local  maximum  in  this  direction.  The  method  is  implemented  as  follows.  Let 

01  =  0^“^^  ;  04  =  0 

02  =  0.618  0  1+0.382  0  4  (2.10) 

03  =  0.382  0  1+0.618  0  4 

and  compute  L  at  these  four  points.  If  there  is  a  single  local  maximum  over  the  segment, 
then  it  can  be  found  by  an  iterative  scheme; 

If  L(0i)  is  the  largest  of  the  four,  replace  04  <—  03  and  return  to  (2.10) 

If  L(04)  is  the  largest  of  the  four,  replace  0i  <—  02  and  return  to  (2.10). 

Failing  these. 

If  L(02)  is  smaller  than  L(0  3),  then  make  the  replacement  0i  <—  02 
If  L(02)  is  larger  than  L(0  3),  then  make  the  replacement  04  <—  03 
and  return  to  (2. 10). 
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Repeat  the  process  until  there  is  no  change.  Then  go  back  to  the  Newton-Raphson 
type  scheme  and  continue  in  a  new  direction.  This  golden  section  augmentation  will  slow 
when  overly  stringent  convergence  criteria  are  used.  But  is  often  preferable  to  the  seeking 
of  a  better  initialization. 

2d.  Likelihood  System;  Complete  Samples 

Let  Xi,  X2,  •••  ,  Xn  be  a  random  sample  of  life  length  random  variables  each  having 
pdf  f(x;  0)  and  survivor  function  S(t;  0).  Note  how  the  parameter  0  is  now  included  in 
the  notation;  it  may  be  multidimensional.  The  likelihood  function  is  expressed 

lik(0)  =  n;i,[f(x.;0)]  and  L  =  x:,log[f(x,;0)].  (2.11) 


The  partial  derivative  of  L,  the  log  likelihood,  with  respect  to  0  has  expected 
value  equal  to  zero  and  that  the  expected  square  of  this  quantity  is  the  negative  of  the 
expected  value  of  the  second  partial  of  the  log  likelihood.  The  requisite  smoothness 
conditions  for  interchanging  the  expectation  operation  with  the  requisite  partial 
derivatives  are  presumed.  The  subscript  0  is  used  to  denote  partial  derivative,  and  the 
format  aspects  of  the  partial  of  the  log  likelihood  take  the  appearance 


Le=Z 


n 

i=l 


f9(Xi;0) 

f(Xi;0) 


(2.12) 


and  the  second  order  partial  derivatives 


T  _V”  f^eeCXjjQ)  rf9(Xii9)-|2^ 
09  ^f(x,;0)^^ 


(2.13) 


where  the  double  subscript  00  refers  to  a  common  component  of  the  vector  0.  Any  mixed 
partial  derivative  of  second  order  has  the  form 


"  I  ^0,02 

^  f(x,;0) 


f9.(Xi;9)f92(Xi;e) 

f^(x.;0) 


}• 


(2.14) 


One  must  develop  these  quantities  for  each  specific  model.  Should  the  system  of 
equations  for  the  mle’s  be  nonlinear,  one  may  use  the  iteration  schemes  described  above. 

The  negative  expectation  of  {Lg  g  }is  nio,  where  lo  is  the  single  observation 
information  matrix.  A  version  of  the  multivariate  central  limit  theorem  says  that 

4n{e-e)  «MLA(0,/o‘) .  (2.15) 
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It  follows  that 


{6  -  6)'  nig  {6-6)~  xlk)  >  where  k  is  the  dimension  of  6.  (2.1 6) 

This  distributional  point  will  be  exploited  in  order  to  find  joint  eonfidenee  ellipses  for  the 
parameters. 

Approximate  joint  eonfidenee  regions  for  the  parameters  ean  also  be  obtained  using 

G"  =  -2[L(^)-L(^')]  (2.17) 


This  too  will  be  used  for  some  eomparisons. 

2e,  Bias  Reduction 


Maximum  likelihood  estimates  are  known  to  be  biased  generally.  It  is  appropriate 
to  reduee  this  bias.  The  teehnique  will  be  presented  using  the  notation  for  eomplete 
samples.  Explieit  relevant  formulae  will  be  generated  where  the  explieit  models  are 
diseussed.  The  use  of  the  teehnique  under  eensored  sampling  follow  these  same  ultimate 
formulae,  but  the  eensored  ease  has  a  mueh  more  diffieult  likelihood  funetion  and  the 
development  of  the  requisite  eomputational  formulae  must  await  another  time.  It  does 
present  a  ripe  area  for  applieation  and  study. 

This  seetion  expands  upon  the  bias  reduetion  analysis  of  mle  estimates  that 
appears  in  [Cox  and  Hinkley,  p.  309].  The  basie  idea  is  to  use  a  third  order  expansion  of 
the  log  likelihood  about  the  mle.  The  teehnique  deals  with  a  single  eomponent  of  the 
parameter  veetor  and  is  expeeted  to  work  best  when  the  estimators  of  those  eomponents 
are  not  strongly  eorrelated. 


The  log  likelihood  funetion  L  is  a  sum  of  n  terms,  where  n  is  the  sample  size. 
Consider  a  single  parameter  0.  The  seore  of  an  observation  is  the  partial  derivative  (of 
one  term)  of  the  log  likelihood  with  respeet  to  0.  The  seore  of  the  i**'  term  (without 
subseript  as  it  plays  no  role)  to  that  sum  will  be  ealled  U;  the  E*  and  2”‘*  partial 
derivatives  will  be  denoted  U'  and  U",  respeetively.  Use  the  dot  subseript  notation  to 
denote  summation  over  the  n  terms.  Thus, 


ain(f) 

50 


f  f  f 

^ee  _  r_!jL]2  ■  pj "  —  see 

f  Lf  J  ,  f 


f  f  f^ 

3^  +  2[^] 


(2.18) 


and  the  expansion  may  be  expressed 


2 


0  =  U.(0)  =  U.(0)  +  (0 -0)U.’(0)  +  -(0  -0)"U."(0  )  +Op(n-‘'0  (2.19) 
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Taking  expectations  through  this  expression,  we  obtain 


E(0-0)E{U.’(0)}  +cov{0,U.'(0)}  +  ^E(0  -0)'E{U."(0)} 

+  ^cov{(0  -  0)^  U."(0)}  =  0(n  (2.20) 

At  this  point  introduce  the  notation,  for  a  single  observation, 

K,^(0)  =  E{[U(0)]^[U'(0)]^},  (2.21) 

and  to  note  that 

E{U"(0)}  =  -3k,, (0)  -  K3^(0)  (2.22) 

and  E{EI."(0)  is  n  times  the  above  amount. 

In  order  to  justify  the  above  expression  and  what  follows  one  should  keep  in  mind  that 
0  =  E(^)  =  E(^)  =  E(^) 

and  cov(0,U'(0))  =  0(l/n)  and  cov[(0  -0)^U"(0)]  =  o(l/n). 

Eurther,  E[U.'(0)  =  ni(0  ),  where  i(0)  is  the  information  scalar  (one  by  one  matrix) 
and  E(0  -  0  )"  «  l/ni(0). 

The  bias  function  is  b(0)  =  E\6  -6  ] .  Now  let’s  take  the  expectation  of  (2.17). 
b(0)  ni(0)  +  1/(2  i(0))  [-3k,  ,  -K3  q]  =  o(l/n) 


and  this  yields  the  approximation 


b(0)  - 


3Kh  +  K30 

2nk(0) 


(This  is  at  variance  with  [Cox  and  Hinkley,  p.  310,  Equation  (35)]). 


(2.23) 
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The  bias  reduction  is  executed  by  using  b{6  ) .  I.e.,  the  reduced  biased  estimator  is 

e*  =  e-b{e).  (2.24) 

2f.  Likelihood  System;  Censored  Samples 

Let  Xi,  X2,  •••  ,  Xn  be  a  random  sample  of  life  length  random  variables  having  pdf 
f(x;0)  and  survivor  function,  S(t  ;0  ). 

Type  I  censoring.  Testing  of  item  i  will  stop  at  f  if  the  item  has  not  failed  by  that  time. 

Type  II  censoring.  All  testing  stops  when  the  r*  item  has  failed,  i.e.,  at  X(r). 

Let  5i  =  1  if  Xi  is  <  (it’s  censoring  time.  Type  I)  f  or  X(r)  (for  Type  II)  (2.25) 

=  0  o.w. 

The  likelihood  function  can  be  expressed  in  a  common  format  for  the  two  kinds  of 
censoring. 


lik(0)  =  nil,  [f  (X. ;  0)]®'  S(t, ;  0)'-«-  (2.26) 

Our  first  goal  is  to  go  through  the  formalities  of  showing  that  the  partial  derivative 
of  the  log  likelihood  has  expected  value  equal  to  zero  and  that  the  expected  square  of  this 
quantity  is  the  negative  of  the  expected  value  of  the  second  partial  of  the  log  likelihood. 
The  requisite  smoothness  conditions  for  interchanging  the  expectation  operation  with  the 
requisite  partial  derivatives  are  presumed.  Let  the  logarithm  of  the  likelihood  function  be 

L(0)  =  ZiljS.  ln(f(x,;0)  +  (l-5,)ln(S(t,;0)],  (2.27) 


and  using  the  subscript  0  to  denote  partial  derivative,  the  format  aspects  of  the  partial  of 
the  log  likelihood  take  the  appearance 


Le 


f9(Xi;0) 

f(x.;0) 


(1-5;) 


SeM)] 

S(t;;0)^ 


(2.28) 


and  the  second  partial  derivatives 


y.  f..(x,;9)  f,(x.;9), 

ee  ^f(x,;0) 


+ 


SeeCljjQ) 

S(t;;0) 


S(t;;0)^  ^ 


(2.29) 
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^  n  fe,e,(x.;9)  _  f9,(Xi;9)fe^(Xi;9) 


f(x.;9) 


fix.;  9) 


"  (1  _ 5  )[  _  Se.(t.;9)Se^(t.;9)^ 


S(t,;9) 


Sit.;  9) 


(2.39) 


The  mle’s  can  be  found  by  setting  (2.28)  equal  to  zero,  and  the  Hessian  can  be  computed 
from  (2.29)  and  (2.39). 


The  expected  values  of  these  quantities  involve  terms  related  to  the  form 
f  (x  ;9)  dx  +  S(t;9)  =  1 , 


(2.31) 


which  when  differentiated  respect  to  9,  produce  the  useful  structures 
J/0(x;9)dx  +  Sg(t;9)  =  9;  J ‘f00(x;9)  =  9 . 


(2.32) 


From  this  we  formally  show  that  the  expected  value  of  (2.28)  is  zero,  and  for 
Equation  (2.29)  we  see  that 


E{L00} = E{  } + z::E(i-5.)[^^] 


f(x.;9) 


S(t.;9) 


and 


EIL  L  1-V"  E18  11^0^192  ,  Y"  nrix 

E{Ee.EgJ  2ji=i  .qm2  1 


[f(x.;9)r 


[S(t.;9)r 


and  E{EgLgJ  =  -E{Ee_eJ. 

The  maximum  likelihood  estimates  will  be  the  same  regardless  of  the  type  of  censoring. 
However,  the  effect  of  censoring  upon  the  information  matrix  and  bias  reduction 
methodology  does  change  with  the  type  of  censoring. 
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3,  Computational  Studies 

This  section  is  likely  the  most  useful  to  the  practitioner.  It  shows  how  to  use  the 
programs  and  suggests  many  kinds  of  exploratory  computations  and  graphs.  The  contents 
are  partitioned  as  follows. 

3a.  Comparison  between  the  Exponential  and  Lognormal  distributions. 

3b.  Comparison  between  the  Weibull  and  Gamma  distributions. 

3c.  Comparison  of  planned  replacement  rules  for  several  models. 

3d.  Methods  for  dealing  with  the  Weibull  distribution;  complete  samples;  bias 
reduction;  censored  samples. 

3e.  Methods  for  dealing  with  the  Gamma  distribution;  complete  samples;  bias 
reduction;  censored  samples. 

3f  Methods  for  dealing  with  the  Lognormal  distribution;  complete  samples;  bias 
reduction;  censored  samples. 

3a,  Comparison  of  the  Exponential  and  Lognormal  Models 

The  exponential  distribution  is  a  favorite  because  of  its  simplicity.  The  lognormal 
distribution  also  appears,  being  derivable  from  fairly  plausible  assumptions  about  certain 
types  of  failure  processes.  [Breiman;  Statistics,  1973,  Houghton-Mifflin,  p.  197]  has 
drawn  attention  to  the  idea  that  these  two  distributions  can  be  adequately  fitted  to  some 
failure  time  data,  and  for  small  samples  it  is  difficult  to  discriminate  between  the  two.  In 
this  example,  the  chi  square  goodness  of  fit  statistic  (D)  is  6.2  with  4  degrees  of  freedom 
for  the  exponential  lit,  estimated  mean  =  41.1.  On  the  other  hand,  the  lognormal 
distribution  produces  a  D  value  of  7.5,  for  p.  =  3.3  and  a  =  1.1.  His  point  is  that  the  chi 
squared  procedure  has  little  discrimination  power,  even  for  a  sample  of  size  50.  In  most 
cases,  the  choice  between  the  two  can  be  resolved  by  comparing  the  two  q-q  plots. 

Ligure  3.1  makes  another  comparison  of  the  two  distributions,  this  time  using 
Exp(l)  and  Lognormal(0,  1).  The  density  functions  appear  similar,  but  the  other  support 
function  allows  one  to  be  more  discriminating  as  to  their  properties.  Generally,  the 
lognormal  hazard  function  has  the  shape  of  an  inverted  “U.”  This  is  implausible  for  most 
situations.  In  spite  of  this  unattractive  feature,  it  has  been  used  in  a  number  of  diverse 
situations.  See  [Lawless,  p.  24]. 

Since  the  hazard  function  for  the  lognormal  is  below  that  of  the  exponential  the 
former  has  a  greater  survivability  and  a  longer  residual  life. 
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Compare  Density  functions  Compare  Survivor  Functions 


0  1  2  3  0  1  2  3 


Compare  Hazard  Functions  Compare  Cum.  Hazard  Functions 


0  12  3 


Figure  3.1 

Behavior  of  the  Lognormal  Model  compared  with  the  Exponential  Model. 

3b,  Comparison  of  the  Weibull  and  Gamma  Distribution  Models 

The  Weibull  distribution  is  a  common  choice  for  a  survival  distribution,  as  is  the 
gamma  distribution.  The  former  is  the  more  popular,  largely  because  the  hazard  function 
follows  a  power  law;  decreasing  for  the  shape  parameter  a  being  less  than  one,  and 
increasing  for  the  parameter  being  more  than  one.  The  hazard  function  for  the  gamma 
law  is  also  DFR  (decreasing  failure  rate)  for  the  shape  parameter  smaller  than  one  and 
increasing  (IFR)  when  it  is  larger  than  one;  but  for  large  values  of  the  variate  it 
approaches  an  asymptote.  Thus,  there  is  a  serious  choice  to  be  made  between  these  two 
laws  even  though  their  density  functions  are  quite  similar.  They  are  both  unimodal  and 
skewed  positively.  It  is  difficult  to  discriminate  between  the  two  based  on  complete 
samples.  The  distinctions  are  in  the  tails.  Some  comparisons  between  the  two  follow. 
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Compare  Density  functions 


Compare  Survivor  Functions 


Compare  Hazard  Functions 


Compare  Cum.  Hazard  Functions 
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Compare  E{Residual  Life} 


Legend 


-  Weibull(0.75,  2) 

mean  =  1.84  stdev  =  2.33 

- Gamma(.72,  .28) 

mean  =  2.57  stdev  =  3.03 


Figure  3.2 

Comparison  of  the  Weibull  and  Gamma  Distribution  Models;  DFR  Case. 

Figure  3.2  makes  the  eomparison  between  the  Weibull  and  gamma  distributions  in 
a  deereasing  failure  rate  ease.  The  Weibull  (  a  =  0.75,  P  =  2)  ease  was  seleeted.  A 
random  sample  of  size  200  was  simulated  and  the  parameters  for  fitting  a  gamma  (a,  X) 
distribution  were  estimated  using  maximum  likelihood.  The  result  is  in  the  legend  of 
Figure  3.2.  The  shape  parameter  a  is  unrelated  to  its  counterpart  for  the  Weibull 
distribution,  but  they  share  the  same  rule  for  discriminating  between  IFR  and  DFR.  The 
parameter  X  is  called  the  rate  parameter. 

Inspection  of  Figure  3.2  shows  that  it  would  be  difficult  to  make  the  choice 
between  these  two  families  without  looking  closely  at  the  tails.  The  effect  of  comparing 
the  expected  residual  life  functions  shows  the  effect  of  a  finite  asymptote  for  the  gamma 
distribution. 
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Compare  Density  functions 
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Compare  Cum.  Hazard  Functions 
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Compare  E{Residual  Life} 
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Legend 


-  Weibull(2,  2) 

mean  =  1.77  stdev  =  2.18 

- Gamma(3.75,  2) 

mean  =  1.88  stdev  =  0.97 


0  1  2  3  4  5 


Figure  3,3 

Comparison  of  the  Weibull  and  Gamma  Distribution  Models;  IFR  Case. 

An  IFR  eomparison  of  these  two  distributions  is  made  in  Figure  3.3.  The  Weibull 
(2,  2)  distribution  was  ehosen  (linear  hazard  funetion)  and  a  reasonable  eonvenient 
matehing  gamma  (3.75,  2)  distribution  was  seleeted  after  some  experimentation.  The 
story  is  mueh  the  same  as  in  the  previous  ease.  The  gamma  hazard  funetion  has  an  ogive 
shape  and  exhibits  a  sharper  separation  from  the  Weibull  hazard  funetion.  The  lower 
hazard  values  translate  in  higher  values  for  the  expeeted  residual  life.  Again,  the  effeet  of 
the  asymptote  is  apparent. 

S-Plus  eodes  for  the  ereation  of  these  graphs  are  in  Appendix  E.  See  the  funetions 

expl .  Inorm  ( )  ,  weil.gamO,  wei2.gam(). 
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3c,  Cost  of  Planned  Replacements 

The  material  in  Seetion  2b  is  used  to  eontinue  our  eomparison  of  the  Weibull  and 
gamma  distributions. 

The  relative  eost  eurves  for  our  two  oases,  Weibull  (2,  2)  and  Gamma  (3.75,  2),  are 
superimposed  for  three  different  ratios  of  o/k,  the  additional  eost  o  of  the  unplanned 
replaoement  to  the  basio  eost  of  a  planned  replaoement.  The  minimum  eost  polioies  are 
tabled; 

Ratio,  o/k  0.5  1.0  2.0 

Min  eost  Weibull(2,2)  0.844  1.091  1.476 

Min  eost  Gamma(  3.75, 2)  0.800  1.043  1.379 

The  members  of  the  paired  eost  eurves  traok  eaoh  other  quite  well.  The  oosts  are 
very  flat  after  the  initial  drop  off.  The  optimum  oosts  are  less  than  one  standard  deviation 
beyond  the  mean.  The  S-Plus  oode  is  oost.oomp(). 

ratio  =  0.5 


3d.  Weibull  Distribution 


Three  types  of  eomputations  are  illustrated:  estimation  for  eomplete  samples;  bias 
reduetion  teehnique  for  eomplete  samples;  and  estimation  eomputations  for  eensored 
samples.  The  model  support  funetions  are  straightforward  and  need  not  be  illustrated. 

Complete  Samples 

The  ball-bearing  data  is  employed,  see  Seetion  1 .  It  is  a  eomplete  sample  of  23 
failure  times  measured  in  millions  of  revolutions.  The  eall 

weibull .  est  (BallB)  (3.1) 

returns  shapeimle  seale:mle  samp  size  shape:mm  sealeimm 
2.101808  81.87422  23  2.01545  81.50336 

and  the  interpretation  is  a  =  shape;  P  =  seale;  mle  =  maximum  likelihood;  and 
mm  =  method  of  moments.  Using  the  mle  values,  the  estimated  information  matrix  is 


J  0.8676723  0.005163827  1 

[0.005163827  0.0006590095J  ' 


(3.2) 


One  may  graph  an  asymptotie  eonfidenee  ellipse  for  (a,  P),  using  the  relationship 

{a-a,p- P)ni{a-a,p- py  ^  x^d)  (3-3) 

for  n  suffieiently  large.  The  funetion  ellipse  (q,  m,  d,  no  =  100)  eomputes  the 

upper  and  lower  portions  of  the  eonfidenee  ellipse.  The  inputs  are  q  =  n/^;  m  is  the 

eentering  veetor  (a,y0)  ;  d  is  the  100(1 -a)**'  quantile  of  the  ehi  square(2)  distribution;  and 
no  is  the  number  of  points  to  use  in  eaeh  quarter  of  the  ellipse.  The  output  is  a  2no  by  3 
matrix.  When  one  plots  the  superimposed  graphs  of  eaeh  of  eolumns  2  and  3  against 
eolurnn  1,  the  result  is  the  ellipse.  We  an  also  plot  an  approximate  eonfidenee  region 
based  on  the  G  distribution.  Equation  (2.17).  The  ealling  sequenee  is 

X  <-  seq(1.4,  2.9,  .1);  y  <-  seq(60,  105,  1) 

zz  <-  Gsq.wei(BallB,  x,  y,  m[l],  m[2],  n)  (3.4) 

plot(plot.wei[,l],  plot.wei[,2],  type="l",  xlab="shape",  ylab="seale", 
xlim  =  e(1.4,  2.9),  ylim  =  e(60,  105)) 
hnes(plot. wei[,  1  ] ,  plot.wei[,3]) 

eontour(x,  y,  zz,  nlevels  =  1  ,v  =  d,  add  =  T,  Ity  =  3  ,labex  =  0) 

title(main  =  "90%  Confidenee  Regions  for  Weibull  Parameters")  (3.5) 

The  results  are  in  Figure  3.5.  The  dashed  eurve  is  eomputed  from  the  G^  distribution. 
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90%  Confidence  Regions  for  Weibull  Parameters 


shape 


Figure  3.5 


Bias  Reduction 


Using  (2.24),  (4.16)  the  formulae  for  bias  reduetion  beeomes 


^  ,  4.309165^  ^  ,  l  +  3/«^ 

a*  =  a(lH - ^ - );p*  =  p(lH — ^ — —) . 


2n 


2nd 


(3.6) 


And  when  applied  to  the  mle’s  found  for  the  ball-bearing  data,  the  reduced  bias  estimates 
become 


a*  =  2.2987,  (mle  =  2.101808) 
p*  =  83.92977  (mle  =  81.87422). 

For  a  sharper  comparison,  let  us  use  a  Monte  Carlo  random  sample  from  a 
specified  Weibull  distribution,  say  Weibull(  1.5,  2).  Let  x  <-  rweibull(15,  1.5,  2).  The 
estimation  function  returns 

weibull.  est(x)  (3.7) 

shapeimle  scaleimle  samp  size  shapeimm  scaleimm 
1.807413  1.720534  15  1.707669  1.707839 
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Use  of  (3.7)  allows  construction  of  the  table 


Parameter 

Actual 

Max,  Lik, 

Bias  Reduced 

a 

1.5 

1.807413 

2.067028 

P 

2.0 

1.720534 

1.804933 

Censored  Samples 

Let  us  treat  the  myeloma  data  set,  per  the  Introduction,  65  observations,  17  of 
which  are  right  censored;  x  <-  myeloma[myeloma$del=  =1,1]  and 
t  <-  myeloma[myeloma$del=  =0,1].  Begin  with  the  finding  of  initial  estimating  values. 
Use  the  call' 

th  <-  mit.wei(x,  t),  (3.8) 

which  returns;  th 

u  b  alph  beta 

3.949889  1.45214  0.6886389  51.92961 

then  the  estimation  function;  use  the  call 

Newt. wei(x,  t,  th[  1:2]),  (3.9) 

which  returns 

u  b  alph  beta  flag 

3.492412  0.9253385  1.080686  32.86512  0 

(Note:  The  flag  starts  at  one  and  is  changed  to  zero  once  the  bisection  search  option  is 
invoked.)  See  Section  4. 

3e,  Gamma  Distribution 

Complete  Samples 

We  use  the  parameterization  (a,  X)  where  X  is  the  data  parameter.  Many  prefer  to 
use  the  scale  parameter,  P  =  1/X.  Because  of  this,  we  carry  the  option  in  the  example. 
Let  us  again  use  the  ball-bearing  data  and  illustrate  the  use  of  the  programs.  The 
immediate  goal  is  to  fit  the  gamma  distribution  and  estimate  the  Information  matrix. 
Model  fit  comparisons  are  made  at  the  close  of  this  section.  Equation  (5.18)  is 
implemented  in  a  function  gam .  es  t  ( )  .  The  application  upon  executing 

th  <-gam.  est  (BallB)  (3.10) 


*  The  parameters  (u,  b)  relate  to  the  Extreme  Value  distribution;  the  bisection  search  method  plays  a  role  as 
well.  These  things  are  explained  in  Section  4. 
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produces  the  return 

alpha-mle  lambda-mle  n  alpha-mm  lambda-mm 
4.024706  0.05572775  23  3.710832  0.05138171 

Note  the  output  contains  both  the  maximum  likelihood  estimator  and  the  method  of 
moments  estimator.  It  follows  that  beta-mle  =  17.94438  and  beta-mm  =  19.46218.  The 
information  matrix  (see  (2.16)),  is  estimated  to  be 


JO.2819  -17.94  1 

[-17.94  1295.959J 


l{a,p) 


[0.2819  0.0557] 

•  (3.11) 

[0.0557  0.0125] 


One  may  graph  an  asymptotic  confidence  ellipse  for,  say  (a,  P),  using  the  relationship 

{a-a,p- p)ni{a-a,p- py  ^ Z(2)  (3-12) 

for  n  suffieiently  large.  The  funetion  ellipse  (q,m,d,  no  =  100)  computes  the 
upper  and  lower  portions  of  the  confidenee  ellipse.  The  inputs  are  q  =  n  / ;  m  is  the 
centering  vector  (a,y0)  ;  d  is  the  100(1 -a)**'  quantile  of  the  ehi  square(2)  distribution;  and 
no  is  the  number  of  points  to  use  in  each  quarter  of  the  ellipse.  The  output  is  a  2no  by  3 
matrix.  When  one  plots  the  superimposed  graphs  of  each  of  columns  2  and  3  against 
column  1,  the  result  is  the  ellipse.  For  an  example  using  the  ball-bearing  data;  n  =  23 
and  d  =  qchisq(0.9,  2);  no  is  100.  Setting  plot.dat  <-  ellipse(q,  m,  d)  and  then 


plot  (plot .  dat  [,  1  ]  ,  plot .  dat  [ ,  2  ]  ,  type  =  "1")  (3.13) 

lines (plot. da t[, 1] ,  plot. da t[, 3] ) 


a  -  X 


a  -  p 


Figure  3.6 

Gamma  Parameters;  90%  Confidence  Ellipses. 
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The  a-X  ellipse  was  eonstrueted  from  the  above  eode.  The  a-p  ellipse  was 
eonstrueted  from  the  former  using  P  =  \/X.  It  is  remarkably  elose  to  the  ellipse  that 
would  be  based  on  the  right-hand  side  of  (3.12). 

Bias  Reduction 

First,  let  us  observe  the  effeet  of  applying  Equations  (5.24),  (5.29),  and  (5.34)  to 
the  ball-bearing  data.  The  only  formal  assistanee  is  offered  in  terms  of  the  S-Plus 
funetion  bias.gam(),  whieh  produees  the  eonstant  C. 

C<-bias.gam(th[l])  =  -0.9936822  (3.14) 

for  use  in  (5.24).  The  results  of  all  three  bias  reduetions  are 

Ball-Bearing  Data 

Parameter  Max,  Lik  Bias  Reduced 

a  4.02471  4.04631 

k  0.05573  0.05633 

p  17.9444  17.9683 

It  may  be  more  useful  to  test  the  method  using  a  smaller  sample  drawn  from  a  Gamma 
population  with  speeified  parameters.  Aeeordingly,  let  us  use 

X  <- rgamma(  15,  shape  =  1.5,  rate  =  0.5).  (3.15) 

The  results  are  in  the  table; 


Parameter 

Actual 

Max.  Lik 

Bias  Reduced 

a 

1.5 

1.2668 

1.2978 

k 

0.5 

0.40196 

0.4231 

P 

2.0 

2.4878 

2.65097 

Censored  Samples 

Let  us  treat  the  myeloma  data  set  (App  E)  65,  observations  17  of  whieh  are  right 
eensored;  x  <-  myeloma[myeloma$del=  =1,1]  and  t  <-  myeloma[myeloma$del=  =0,1]. 

Eirst,  we  need  initial  estimates  for  input  into  the  iteration  method.  The  funetion 
init.gam(x)  takes  the  uneensored  portion  of  the  data  and  eomputes  the  method  of 
moments  estimates.  E.g., 

thO  <- mit.gam(x)  and  the  return  is  0.98106187  0.04014575.  (3.16) 

Next,  we  use  the  general  funetion  Itest(x,  t,  thO,  Newt,  Llik),  whieh  will  exeeute 
the  iteration  method  deseribed  in  Seetion  2o.  The  inputs  are 
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X  the  uncensored  data 
t  the  censored  data 
thO  the  initial  estimates 
Newt  <-  Newt.gam 
Llik  <-  Llik. gam 

The  entry  Itest(x,  t,  thO,  Newt.gam,  Llik.gam)  (3.17) 

returned 


d  =  1.06707553;  ^  =  0.03315037;  number  of  cycles  to  convergence  =  12. 
The  codes  for  this  and  support  functions  can  be  found  in  Appendices  C  and  D. 

3f.  Lognormal  Distribution 
Complete  Samples 

Let  us  fit  the  lognormal  distribution  to  the  ball-bearing  data.  For  complete 
samples  this  is  a  very  simple  task.  Set  x  <-  BallB  and  then 

ju  <- mean(log(x))  =  4.150383;  d  <- sqrt((n-l)/n)*stdev(log(x))  =  0.5216865.  (3.18) 


Let  us  also  generate  the  90%  approximate  confidence  ellipse  for  (p,  a),  see 
Equation  (2.16)  and  compare  it  with  other  confidence  regions  for  normal  data, 
specifically  a  region  based  on  of  Equation  (2.17),  and  an  exact  trapezoidal-shaped 
region.  These  are  explained  in  Section  6b.  See  Equations  (6.16)  and  (6.17). 

Eor  the  former  we  require  the  estimated  information  matrix 


J3.674352 

1  0 


“  !■ 

7.348704] 


(3.19) 


and  then  m  <-  c(4. 150383,  0.5216865);  d  <-  qchisq(0.9,  2);  plot.dat  <-  ellipse(n*I,  m,  d). 
The  plotting  sequence  for  Eigure  3.7  is 

plot(plot.dat[,l],  plot.dat[,2],  type="l",  xlab="mu",  ylab="sigma",  ylim=c(.3,.8), 
xhm=c(3.8,  4.5)) 
hnes(plot.dat[,l],  plot.dat[,3]) 

contour(x,  y,  G2,  nlevels  =  l,v  =  d,  labex  =  0,  lty=3,  add=T) 

(see  Section  6.b  for  the  computation  of  G2) 
trap  <-  NormCT(log(BallB),  0.1,  graph=E) 
hnes(trap[l,],  trap[2,]) 
points(xb,  s) 

title(main  =  "90%  Confidence  Regions  for  Mu  and  Sigma") 
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90%  Confidence  Regions  for  Mu  and  Sigma 


mu 


Figure  3,7 

The  solid  ellipse  is  based  on  (2.16).  The  dashed  curve  is  based  on  the  G  statistic 
of  (2.17).  The  trapezoid  region  is  an  exact  region  based  on  the  independence  of  the 
sample  mean  and  sample  variance. 

Bias  Correction 

The  estimator  for  p.  is  unbiased  and  the  estimator  for  a  need  only  be  divided  by 
the  correction  factor  in  (6.18).  For  the  ball-bearing  data  this  results  in 

[i*  =  4.150383  and  a*  =  1.034156x0.5216865  =  0.5395052  (3.20) 

Censored  Samples 

Use  the  myeloma  data  from  Appendix  E.  The  initializing  values  are  taken  from 
the  mean  and  standard  deviation  of  the  log  of  the  uncensored  data.  Call  it  thO.  I.e., 

thO  <-_c(mean(log(x)),sqrt((22/23)*var(log(x)))),  (3.21) 

which  returns  4.1503827  0.5216865. 

Then  execute 

Paramest  <-  Itest(x,  t,  thO,  Newt.logn,  Llik.logn)  (3.22) 
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and  the  return  is 

ju  =  4.2533626;  a  =  0.5216074;  number  of  eyeles  =  18. 

Model  Selection 

Having  maximum  likelihood  estimates  for  the  parameters  of  our  eompetitive 
models  for  deseribing  the  ball-bearing  data,  it  is  interesting  to  eonsider  how  one  might 
ehoose.  In  [Meeker  and  Eseobar]  the  ehoiee  is  an  extended  gamma  distribution,  a 
distribution  that  is  not  direetly  treated  here.  It  is  easy  to  eompute  the 
Kolmogorov-Smirnov  test  statisties  and  use  them  as  distanee  funetions.  Meaningful 
p-values  eannot  be  eomputed  beeause  the  distributions  are  fitted  from  data.  The  distanees 
are 


ks.distanee(BallB,  Weibull)  =  0.151 
ks.distanee(BallB,  Gamma)  =  0.123 
ks.distanee(BallB,  Lognorm)  =  0.090 

The  ball  bearing  set  is  a  eomplete  sample  and  we  ean  make  a  graph  that  eontains 
the  empirieal  and  the  fitted  model  distributions.  For  an  empirieal  distribution  we  use 

=  jlin  +  \)  for  j  =  1,  -  ,  n  (3.23) 

and  plot  these  diserete  values  against  the  order  statisties  . 


CDF's  of  the  Data  and  the  Three  Fitted  Models 


Figure  3,8 
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4,  The  Weibull  Distribution 


The  Weibull  distribution  is  a  very  popular  model  for  reliability  work,  largely 
beeause  of  its  ease  of  use  and  of  its  monotone  hazard  funetion,  whieh  is  inereasing  if  the 
shape  parameter  a  is  larger  than  one  and  deereasing  when  that  parameter  is  smaller  than 
one.  The  parameter  P  is  the  seale  parameter.  The  life  length  random  variable  X  has  a 
Weibull  distribution  if  Y  =  (X/p)“  has  an  Exp(l)  distribution.  This  relationship  is 
exploited  broadly  in  what  follows. 

The  five  seetions  in  this  ehapter  deseribe  formulae  for:  a)  the  model 
eharaeterization  funetions;  b)  likelihood  analysis  for  eomplete  samples;  e)  bias  reduetion 
of  maximum  likelihood  estimates;  and  d)  likelihood  analysis  for  eensored  samples.  This 
last  section  includes  details  of  the  use  of  the  extreme  value  distribution  and  its  role  in  the 
estimation  problems. 

4a,  Model  Characterization  Functions 

The  density  and  survivor  functions  are 

f(x)  -  |(^r'expH|n  and  S(t)  =  expHV]  for  x,>0,  t>0,  (4.1) 

where  a  >  0  is  a  shape  parameter  and  P  >  0  is  a  scale  parameter.  The  mean  and  variance 
are 


^  =  p  r(l  +  1/a)  =  p^  [r(l  +  2/a)  -  r^(l  +  1/a)].  (4.2) 

The  hazard  function 

h(t)  =  -(-)“-'  .  (4.3) 

p  p 


The  cumulative  hazard  function 

H(t)=(^)“.  (4.4) 

The  integrated  survivor  function 

SS(t)  =  £- j  "  y''“V%  =  £-{r(l/a)-IncGam[(t)“].  (4.5) 

a  •'(p)  a  P 

The  mean  residual  life 

LL(t)  =  SS(t)/S(t).  (4.6) 
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4b,  Likelihood  Analysis  for  Complete  Samples 

Let  Xi,  X2,  •••  ,  Xn  be  a  eomplete  random  sample  from  a  parent  population, 
X  ~  Weibull(a,  P).  Let  i  be  the  log  likelihood  of  a  single  observation,  x.  The  direet 

analysis  of  this  likelihood  system  follows.  Liberal  use  is  made  of  the  faet  that  =  T  ~ 
Exp(l). 


I  =  ln(a)  -  aln(P)  +  (a-l)ln(x)  -  (— )“ 


(4.7) 


=  --ln(p)  +ln(x)  -  (J)“ln(J  )  =  i[l  +  ln(T)-nn(7)] 
a  P  P  a 


=  -[  -1+  (-)“]  =  -[-1  +  L] 

P  pL  Vp2 


(4.8) 


1 


1 


faa  =  -  ^[  1  +  (-)“ln^(T7)“]  =  -^[l  +  Lln^(L)] 


a 


P  P 


a 


pp  =  l-(a  +1)(|)“]=  ^[l-(«  +  l)T] 


(4.9) 


fap  =  -  [  -!  +  (-)“ +  (-)“ln(-)“]  =  -[-l  +  7  +  nn(7)] 

P  p  L  Vp^  Vp^  Vp^  p 

Sinee  0  =  E{f„}  =  E{  f  p},  we  may  deduee  some  interesting  relationships. 


E{  (jT}  =  E{Y}  =  landE{(^)“ln(^)}-E{ln(X)}  =  --ln(p)  (4.10) 


a 


Next,  when  we  replaee  x  with  Xi  and  sum  over  the  data,  we  may  write  the 
members  of  the  Hessian  as 


La  a 


L„=  ^{n  -  (a  +1)L:.  (f  n 


.X;  , 


Up  =  (yr+z:  (^)“in[(|-r]} 
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The  maximum  likelihood  estimates  are  eomputed  using  an  iteration  funetion  that 
follows  readily  from  the  pair  of  equations  La  =  0  and  Lp  =0.  These  two  equations  are 


easily  extraeted  from  (4.9).  The  latter  yields  the  equation 


[  (~T  ]  1 5  whieh,  when 

b 


substituted  into  the  former,  produees 


1  _1 
a  n 


(4.11) 


This  in  turn  allows  a  determination  of  a  from  the  left-hand  side  after  an  initial  value  of  a 
is  plaeed  into  the  right-hand  side.  The  iteration  proeeeds  when  the  new  value  so 
determined  is  inserted  into  the  right-hand  side  and  the  proeess  repeated.  This  is  often 
ealled  the  natural  iteration  funetion;  it  need  not  eonverge  in  general,  but  in  this  ease  it 
does.  Its  use  is  illustrated  in  Chapter  3. 

The  information  matrix  is 


j  ^  |[l  +  'F'(2)  +  'L"(2)]/a'  ^{2)1/3} 

^(2)1  p  [aipf]' 

Proof  It  follows  from,  see  Appendix  A, 

E(j)“ln"[(A“]=  EjYInW)!  =  r"(2) 

E(A“ln[(A“]  =  EjYln(Y)}  =  r(2) 

r'(2)  =  ^^(2)  and  r"(2)  =  (2)  +  [T^(2)]l  □ 

Use  of  this  is  made  in  (3.3). 

4c.  Bias  Reduction 

The  teehnique  being  used  treats  the  two  parameters  separately, 
(i)  shape  parameter,  0  =  a,  so  U  =  f  a  when  referring  to  (2.18). 

U  =  -{l  +  ln(T)-nn(T)} 
a 


U'  =  \[\  +  Y\n\Y)-\ 
a 
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First, 


^  1 

E{u^}  =  —  {[1 +  3r(i)-3r(2)]  +3[r"(i)-2r"(2)  +  r"(3)] 
a 

+  [  r"’(l)  -3r"’(2)  +  3r’"(3)  -  r"'(4)]}  and  the  expectation  is 

K30  =  ^  {-7.921007}  (see  Appendix  A).  (4.13) 

a 

Second, 

E{uu'}  =  ^{i+r'(i)-r'(2)+r"(2)+r"'(2)-r"'(3)} 
a 

and  the  expectation  is 

Kii  =  ^{-2.136823}  (see  Appendix  A).  (4.14) 

a 

Third,  from  (4.1 1)  the  information  in  an  observation  is 

i(a)  =  lo[l,l]  =  [1  +  ‘F'(2)  +  'F"(2)]/«'  =  ^{1.823681}.  (4.15) 

a 

It  follows  from  (2.20)  that  the  bias  function  is 

b(a)  =  —(-4.309165).  (4.16) 

In 

(ii)  scale  parameter,  0  =  PandU  =  f  p  when  referring  to  (2.15). 

U  =  -[Y-1} 

P 

U'  =  ^[l-(«  +  l)7] 

K30  =  E{U^}  =  (^)'E{Y'-3Y' +3Y-1}  =  2(^)'  (4.17) 

K„  =  E{UU'}  =  ^E[\Y-\-\\\-{a  +  m)  =  (4-18) 
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(4.19) 


i(p)  =  Io[,2,  2]  =  (a/p)' 

It  follows  from  (2.20)  that  the  bias  function  is 

b(p)  =  ^(l  +  3/«).  (4.20) 

2na 

4d.  Treatment  of  Censored  Data;  Use  of  the  Extreme  Value  Distribution 

The  programs  offered  for  complete  samples  use  the  Newton-Raphson  iteration 
scheme  in  two  dimensions.  However,  if  one  transforms  the  data  by  the  logarithm  the 
resulting  density  function  is  that  of  the  extreme  value  distribution.  The  advantage  of  so 
doing  allows  the  elimination  of  one  of  the  two  parameters  in  the  system  of  likelihood 
equations.  This  technique  is  employed  in  treating  the  censored  case. 

The  development  of  the  relationship  between  the  Weibull  and  Extreme  Value 
distributions  is  essentially  that  appearing  in  [Lawless].  Begin  with  the  survivor  function 

S(t;a,  p)  =  Pr{X>t}  =  exp  {  -  (t  /  p)“} .  (4.21) 

It  follows  that  the  pdf  is 

(4.22) 

Both  parameters  are  positive;  a  is  the  shape  parameter  and  P  is  the  scale  parameter. 

The  pdf  of  the  extreme  value  distribution  has  pdf 

1  (— )  (— ) 

*  exp{-e  *  }  for  -go  <  v  <  go.  (4.23) 

b 

This  distribution  is  related  to  the  Weibull  by  the  transformation 

V  =  log(X)  X  =  exp(V)  (4.24) 

and  the  parametric  identification 

u  =  log(p)  p  =  exp(u) 

b  =  1/a  a  =  1/b,  (4.25) 

V  —  X 

from  which  it  follows  that -  =«log( — )  =  log[(x/p)“]. 

b  (5 
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Upon  finding  the  mle’s  for  the  extreme  value  distribution  and  using  the  invarianee 
property  of  maximum  likelihood  estimators  one  can  obtain  maximum  likelihood 
estimators  for  the  Weibull  distribution.  The  development  will  not  depend  on  whether  the 
censoring  mechanism  is  of  Type  1,  Type  2,  or  any  other  right  censoring  plan. 

Let  Xi,  X2,  •••,  Xn  be  the  life  lengths  of  n  items  placed  on  test.  Under  Type  2 
censoring,  testing  stops  when  r  items  have  failed.  Of  course,  if  r  =  n,  then  there  is  no 
censoring.  Under  Type  1  censoring  each  expiring  has  a  test  termination  time  t.  The 
structure  of  the  likelihood  system  allows  both  cases  to  be  treated  with  a  single  set  of 
equations.  The  set  D  contains  those  life  lengths  that  were  completed  prior  to  the 
termination  time.  The  set  C  contains  those  life  lengths  that  exceeded  the  allotted  time  The 
structure  of  the  likelihood  equations  is  given  in  Equations  (2.27)  and  (2.28).  The  indicator 
variables  5i  tells  us  that  Vi  =  In(xi)  when  equal  to  one,  and  Vi  =  In(ti)  when  equal  to  zero. 
The  development  should  be  compared  with  [Lawless,  eq.  (4.1.1)].  The  log  likelihood 
function  has  the  structure 

L(u,b)=  y„ln(/(v,))+X,lii(S,(v,)),  (4.26) 

where  f  has  the  form  (4.22)  and 


V  —  1J 

Sv(u)  =  Pr{V,>v,}  =  Pr{X>h}  =  exp[-exp(^)]  (4.27) 

b 


L  =  ln//k(M,h)  =  -rln(h)  +  ^^^^^^^ 


w.-  u . 


where  ^  exp[^^]  = 


4= 

b  b  I  b 


^  r  1  V-l  V.-U  I'^.V;  -M  W;  -M, 

- >  - +  -  >  - exp(-^ - )  . 

b  b^^  b  b^  b  b 

Setting  the  first  of  these  two  equal  to  zero  and  solving  leads  to  the  separation 


e-=  [tyexp(U]* 

r  b 


(4.28) 


and  substituting  into  the  second  equation  leads  to  the  nonlinear  equation  in  b. 
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(4.29) 


h(b)  =  X'>’,exp(^)/2‘expA)-*-i^„v,=0. 
b  b  r 


This  can  be  solved  using  Newton-Raphson.^  To  this  end,  we  reeord  the  derivative 


Z*  2  vJb  rV^*  v,-/6n2 

iiVu; - ^ +  - ^-1. 


(4.30) 


Having  b  one  obtains  u  from  (4.28);  one  eonverts  to  (a,  P)  using  (4.26). 

Equation  (4.27)  ean  provide  a  way  to  get  initial  estimates  for  the  parameters;  it  works  for 
the  eomplete  ease  as  well.  Speeifieally,  use  the  Kaplan-Meier  estimate  for  the  survivor 
funetion  of  V,  the  extreme  value  variate.  Note  that 

y  =  ln(-ln(5(v))  =  (v -u)/b.  (4.31) 

The  least  square  estimates  of  u  and  b  ean  serve  for  initialization. 


^  We  have  experienced  lengthy  oscillation  using  this  method.  Success  has  been  achieved  by  using 
Newton-Raphson  until  two  consecutive  values  of  h  are  of  opposite  signs.  At  this  point,  we  switch  to  a 
bisection  search.  The  output  contains  a  flag,  which,  if  equal  to  zero,  tells  us  that  the  bisection  search  was 
invoked.  See  the  function  Newt.wei(). 
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5.  Gamma  Distribution 


The  gamma  distribution  mimics  the  Weibull  distribution  in  the  central  portion,  but 
there  are  major  differences  in  the  tails.  The  gamma  distribution,  for  shape  parameter  a 
more  than  one,  has  an  S-shaped  hazard  function  which,  although  monotone  increasing, 
approaches  a  finite  asymptote.  In  this  development,  the  rate  parameter  is  X. 

5a,  Model  Characterization  Functions 

The  density  and  survivor  functions  are 

f(x)  =  — — x“  'exp[-A,x]  and  S(t)  =  l-IncGam(A,  t)  for  x,>0,  t>0,  (5.1) 

r(a) 

where  a  >  0  is  a  shape  parameter  and  A,  >  0  is  a  rate  parameter.  The  mean  and  variance 
are 


|u  =  a/X  =  a/X^.  (5.2) 

The  hazard  function  is 

h(t)  =  f(t)/S(t).  (5.3) 

The  cumulative  hazard  function  is 

H(t)  =  -log[S(t)].  (5.4) 

The  integrated  survivor  function  is 

SS(t)  =  Siu)du  =  uS{u)\"  udSiu) 

a  -u 

=  -tS(t)  +  f  uf  (u)  du  =  - 1  S(t)  +  a  f  — -  du 

•'t  r(a  +  l) 

=  -tS(t,  a)  +  aS(t,  a+1).  (5.5) 

The  mean  residual  life  is 

LL(t)  =  SS(t)/S(t)  (5.6) 

5b,  Likelihood  Analysis  for  Complete  Samples 

Let  Xi,  Xi,  •••  ,  Xn  be  a  complete  random  sample  from  a  parent  population, 
X  ~  Gamma(a,  X).  The  log  likelihood,  scores  and  their  partial  derivatives  have  the  forms 
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I  =  a  ln(A,)  -  ln(r(a)  +  (a-l)ln(x)  -  A,x  (5.7) 

=  ln(}i)  - 'F(a)  +  ln(x)  (5.8) 

I =  a/X  -  X  (5.9) 

=  -na);  =  -a/^^  £„,  =  1/X  (5.10) 

But  if  the  alternative  parameterization  is  used,  P  =  1/A,,  one  must  adjust  (5.9)  and  (5.10) 
and  use 

=  -a  /  P  +  X  /  P^;  £^^  =  -1  /  P^;  £^^  =  a  /  P^  -  2  x  /  P^  (5.11) 

Since  0  =  E{  £^  }  we  may  deduce  an  interesting  relationship. 

E{X-ln(X)}  =  ln(}.)  - 'F(a)  (5.12) 

It  is  useful  to  record  the  statistics 

X  =  Exi/n;  ln(x)  =  E  ln(xi)/n;  s^  =  E(xi  -  x)^/(n  -  1),  (5.13) 

and  then  the  log  likelihood  function  is  easily  expressed  as 

L  =  n  •  a  •  ln(A,)  -  n  •  ln[r(a)]  +  (a  -  1)  E  In(xi)  -  A,  E  Xi,  (5.14) 

and  the  two  components  of  the  gradient  vector  are 

La  =  n  [  ln(A.)-r'(a)/r(a)] +E  In(xi)  (5.15) 

Lx  =n[a/A,]-Exi  (5.16) 

Setting  La,  =  0  allows  the  substitution  of  a  =  Xx  into  La  =  0.  The  resulting  equation 
can  be  solved  by  Newton-Raphson  iteration. 

The  Newton-Raphson  algorithm  requires  initialization  estimates.  It  is  convenient 
to  use  the  method  of  moments  estimators  for  this  purpose.  Thus,  we  solve  the  equations 

X  =  a/X  and  s^  =  a/X^,  (5.17) 

and  obtain  a  =  x  /s  and  A,  =  x  /  s  . 

Because  of  the  elimination  of  A.,  one  requires  only  a  for  initialization  into 

g(a)  =  \)/(a)  -  ln(a)  -  ln(x)  +ln(x)  =  0  (5.18) 

g'(a)  =  -  1  /  a. 
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This  is  managed  by  the  gamma  estimation  function  suite  (see  gam. ml e  ( )  .  Its  use  is 
illustrated  in  3e).  The  S-Plus  code  is  in  Appendix  E.  Asymptotic  expansions  for  the  psi 
function  and  its  derivatives  are  recorded  in  Appendix  A. 


The  information  matrix,  developed  in  Appendix  C,  one  for  each  parameterization 

Up  1 


l{a,X)  = 


\y/'{a)  -\IX\ 

[-yx  ajX}] 


I{a,p)  = 


Up  a!  p^ 


(5.19) 


These  developments  are  applied  in  3e  to  real  data.  The  computations  include 
confidence  regions  for  the  two  parameterizations. 


5c.  Quantities  Needed  to  Execute  the  Bias  Reduction  Method 


The  general  formulae  (2.17),  (2.20),  and  (2.23)  are  developed  for  when  dealing 
with  complete  samples  from  the  gamma  population.  The  method  adopted  treats  the 
parameters  individually.  Accordingly,  the  shape  parameter  a  is  treated  first,  and  then  the 
rate  parameter.  A,. 


Case  i.  Shape  parameter;  0  =  a  and  hence  U  =  f  „ 


(5.20) 

(5.21) 

(5.22) 

-5 

'E(a)}  ,  where  Y  ~  gamma(a,  1) 

(5.23) 

E{U'U}  =  0  since  U' is  constant. 

It  follows  that 


U  =  ln(}t)  - 'E(a)  +  ln(X) 

U'  =  -'E'(a) 
i(a)  =  'E'(a) 

E{U^}  =  E{ln(}.X)-'E(a)}^  =  E{ln(Y)- 
=  t^"(«) 


h{a)  = 


T^"(«) 

2n[^>'{a)f 


(5.24) 


and  the  bias  reduced  estimate  is 


a=a-Clln  with  C  =  ^  \  (5.25) 

['E'(«)]^ 

Accordingly,  the  quantity  C  must  be  computed  case  by  case.  The  S-Plus  function 
bias.gam()  in  Appendix  C  can  be  used  to  compute  it. 

Case  ii.  Rate  parameter;  0  =  A,  and  hence  U  =  f  a. 


34 


U  =  --X,U'  =  =  alX^. 

E{U^}  =  -2alX^ 

E{UU’}  =  0. 

It  follows  that 
b(X)  = 

na 

and  the  bias  reduced  estimate  is 

r=i(i+^) 

nd 

Case  iii.  (Alternate  parameterization),  0  =  A,  and  hence  U  =  f  p 
U  =  X/p^-a/p 
U'  =  -(2X-ap)/p^ 
i(p)  =  (a/p)'. 

Using  Y  =  X/p  leads  to  E{Y-a}^  =  2a  and  E{(  Y-a)(  2Y-a)}  =  2a 
K30  =  2a/p^;  Kij  =  -2a/p^  and 

b(p)  =  -  2  p  /  na^ 

and  the  bias  reduced  estimate  is 

p*  =  y0(l+2/n«'). 


(5.26) 

(5.27) 

(5.28) 


(5.29) 


(5.30) 


(5.31) 

(5.32) 

(5.33) 


(5.34) 

(5.35) 


5d,  Treatment  of  Censored  Samples 

Equations  (2.25)  thru  (2.27)  take  the  following  forms  when  sampling  from  a 
censored  gamma  population. 

E  =  r  [alog(k)-v(a)] +  (a-l) 


+  Xclog[S(ti)] 


(5.36) 


Ea  =  r[log(k)  -\|/(a)]+  ^^log(Xi)  +  ^^S„(ti)/S(ti) 
E,  =  ra/k-X^x,  +  XcS.(t.)/S(h) 


(5.37) 
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(5.38) 


L„„  =  -rv|/'(a)+  Xc(S..(t.)/S(ti)-[S.(ti)/S(ti)f} 

Ui  =  -va/X-+  X;c<S>.i(‘i)'S(ti)-[Si(ti)/S(t,)f} 

L„i  =  r/A+  Xc(S.i(t|)/S(t.)-S„(t,)S,(t,)/[S(t,)f} 

The  unresolved  computational  problems  faced  when  dealing  with  Equations 
(5.37)  and  (5.38)  appear  in  those  summations  over  the  set  C.  They  contain  derivatives  of 
the  Incomplete  Gamma  function  and  pose  a  major  development.  Our  first  step  is  to 
re-express  these  quantities  in  terms  of  he  survivor  function  of  the  standard  gamma 
random  variable  (A,  =  1).  Let  us  use  the  notation 

S*{t)  =  — — f  x“  'e  "'d'x.  (5.39) 

It  is  easily  seen  that  the  above  partial  derivatives  are  needed  only  for  this  standard  form. 
I.e., 


S{t)  =S\Xt)  ■  S„{t)  =  SliXt)  ■  S^{t)  =  SliXt)  (5.40) 

and  the  algorithms  used  are  described  in  Appendix  C.  But,  for  the  mixed  partial 
derivative  and  the  ones  with  respect  to  X,  we  develop  the  formulae  contained  in  the 
summary  below. 

Density  (5.41) 

P(x)  =  X  e  /  r(a)  is  used  for  the  standard  form,  i.e.,  A,  =  1. 
f(x;  a.  A,)  =  — — x“  '  e  =  A,  f  *(A,x)  (or  f(x)  for  short) 

r(«) 

fKx)  =  f(x)[— -1] 

X 

f «  (x)  =  f(x)  [ln(}tx)  -  'F(a)] 
fi,(x)  =  (a/A.  -  x)  f(x) 

f«i(x)  =  [(a/A,  -  x)  (ln(A,x)  -'F(a))  +1/A,]  f(x). 


Survivor 


(5.42) 
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/•OO  *  *  /"OO  * 

S(t;  a,  A,)  =  I  f  (x)dx  =  S  (A,t),  where  S  (t)  =  J  f  (x)dx 
Sl(t)  =  f’(x)[ln(x) -'F(a)]dx 

S„(t)  =  S:(Xt) 

S,(0  =  -^l^(^t)  =  -4f(t) 

A 

S„.(t)  =  -tfUt)[ln(^t)-'F(a)]  =  -4f(t)[ln(Xt)-'F(a)]. 

A 

Systems  of  seeond  partial  derivatives 

f_(x)  =  f(x){[log(Xx)-'P(a)f-na)} 

faA  (x)  =  f(x){(a/?t  -x)  [log(?.x)  -'F(a)]  +  1/?^} 
fiA(x)  =  f(x){[a/?i-xf 

s„„(t)  =  sL(xt) 

SL(t)  =  f  r(x)[ln(x)-'F(a)fdx-'F'(a)S(t) 

S„.(t)  =  -  tf  (^t)[ln(Xt)-'F(a)] 

S,,(t)  =  -t2f(^t)[^-l] 

At 


(5.43) 
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6,  The  Lognormal  Distribution 

The  lifetime  X  is  lognormal,  LN(|u,  a^)  if  Y  =  log(X)  is  N((|a.,  a^).  The  model 
can  be  derived  from  fairly  plausible  assumptions  and  often  found  suitable  for 
representing  lifetimes,  especially  when  large  values  are  not  of  interest.  Some  applications 
are  cited  in  {Lawless,  p.  24],  This  distribution  has  some  strange  properties;  e.g.,  its 
hazard  function  is  an  inverted  bathtub. 


There  is  great  support  for  the  normal  distribution,  and  this  support  translates  to  the 
lognormal  distribution.  For  purposes  of  leveraging  this  support  it  is  convenient  to  draw 
attention  to  some  properties. 

Properties  of  the  Normal  Distribution 

Use  (p  (x)  for  the  N(0,  1)  density  and  0(x)  for  its  cumulative  distribution  function. 

That  is 


(p  (x)  =  exp(-^) 

V2;z-  2. 


for  the  standard  normal  pdf  and  the  properties: 


(6.1) 


r  ^ 

(p'  (x)  =  -  X  (p(x);  0(x)  =  (p(v)dv. 

J  -00 


(6.2) 


Further,  the  survivor  function  for  the  N(0,  1)  case 
S(t)  =  l-O(t); 

and  the  integrated  survivor  function  can  be  developed 


(6.3) 


/•OO  /•OO  /"OO  /•OO  CV 

SS(t;  0,  1)  =  [  [l-0(x)]dx  =  [  [  (p  (v)dv  dx  =  [  (p(v  )[  [  r/xjdv 

Jt  Jt  Jx  Jx  Jt 

/•OO  /"OO  /"OO 

=  J  (v-t)(p(v)dv  =  “J  (p'(v)dv-tj  (p(v)  dv 


=  (p(t)-t[l-0(t)]. 

The  general  N(p,  a^)  case  can  be  expressed  using 

f(x)  =  —  (p(^^)  and  S(t)  =  (p(v)Jv 

a  a 


(6.4) 


(6.5) 
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and  the  integrated  survivor  funetion  is 

SS(t)  =  o  {<p  (^)  -  (^)[  1  -  (6.6) 

a  a  a 

So  mueh  for  the  normal  distribution.  Let  us  turn  to  the  lognormal  distribution. 

6a,  Model  Characterization  Functions 

The  density  and  survivor  funetions  for  LN(|u,  a  )  are 

f(x)  =  —  (p  and  S(t)  =  1  -  0(  )  for  x,  >0,  t>0,  (6.7) 

<JX  <J  (J 

where  |u  and  a  >  0  take  their  usual  interpretation.  The  mean  and  varianee  are 

mean(X)  =  exp[//  +  (T^ /2]  var(X)  =  exp[2p+G^]  [exp(a^)  -  1].  (6.8) 

The  hazard  funetion,  or  age  speeifie  failure  rate  is 

h(t)  =  f(t)/S(t).  (6.9) 

The  eumulative  hazard  funetion  is 

H(t)  =  -log[S(t)].  (6.10) 

The  integrated  survivor  funetion  is 

SS(t)  =  [1  -t}.  (6.11) 

<j 

The  mean  residual  life  is 

LL(t)  =  SS(t)/S(t).  (6.12) 

Only  the  integrated  survivor  funetion  is  a  bit  intrieate  to  understand.  Let’s  first  treat  a 
related  integral. 

/•OO  1  1*00 

[  [1  -  0(w)]  e'^  dw  =  -  [  [1  -  0(w)]  d(e'^"') 

i  a  CT  ^  ^ 

I  <*00 

=  -  {  -  e'^"  [1  -  0(a)]  +  [  e'"Xw)dw} 

rr 
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and  use  the  eomplete  the  square  technique  in  the  integral  part 


1  2 
—  {  -  e'^“  [1  -  0(a)]  +  [  (p(w-G)dw}. 

cr 


With  this  in  mind,  let  us  turn  to 

SS(t)  =  r  [1  -0(  =  ae^  f  [1  -  ®(w)]  e'""' dw, 

Jt  Q-  Ja 


where  a  =  (log(t)  -  p)/a.  Now  we  put  in  the  related  part  and  complete  the  function  in  a 
computable  form. 

SS(t)  =  e^  {  -  o'""  [1  -  0(a)]  +  f  ”  (p(w)dw} 

J  a-a 


[1  -  O  (a  -  0-)  ]  -  t  [1  -  0(a)]. 


□ 


(6.13) 


6b,  Complete  Samples 

Let  Xi,  X2,  •••,  Xn  be  a  complete  random  sample  from  a  parent  population, 
X  ~  LN(|u,  a^).  The  relationship  connecting  Y  =  log(X)  does  not  involve  the  parameters. 
So,  the  well-known  formulae  for  the  maximum  likelihood  estimates  of  p,  and  a  using  yi, 
y2,  •••  ,  yn,  the  logarithms  of  the  data,  provide  the  maximum  likelihood  estimates.  The 
invariance  principle  is  invoked.  Further,  the  information  matrix  will  be  the  same. 


Io(p,  a) 


'2 

cr" 

0 


(6.14) 


Note:  Had  we  been  estimating  instead  of  a,  the  lower  right  entry  would  have  been 
1/(2g^). 

Since  the  normal  and  log  normal  distributions  have  the  same  set  of  parameters, 
then  the  joint  confidence  regions  are  the  same.  The  approximate  ellipses  can  be  generated 
in  the  same  way  as  before;  but  an  exact  trapezoidal-shaped  confidence  region  is  also 
available,  and  it  is  interesting  to  compare  the  two.  The  former  can  be  found  by  using 
(6.14)  in  Equation  (2.16).  The  latter  is  based  on  the  independence  of  the  sample  mean  and 
variance.  Let 


Y  =-yY.  and  SS  =  yiY.-Yf 
n 

when  Yi,  Y2,  •  •  • ,  Yn  is  a  random  sample  from  N(p,  a  ). 


(6.15) 
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The  former  is  distributed  aeeording  to  N(|u,  o^ln)  and  the  latter  by  cr^ xln-\)  ■  It  follows 
that  a  joint  eonfidenee  region  for  (|u,  a)  ean  be  obtained  from  setting  the  produet 


Vx{-z^<4n{^^^)<z^)Vx{k,<SS !  cj^  <k^)  =  1  -  a.  (6.16) 

<j 

Let  the  first  faetor  have  probability  1-ai  and  the  seeond  l-ai.  A  default  position  might 
be  to  use  (l-a)"^  for  eaeh,  but  it  is  useful  to  be  more  flexible.  We  ehoose  to  use  a 
parameter  p,  0  <  p  <  1,  sueh  that 

1-ai  =  (1-a)^  and  l-ai  =  (1-a)^^^. 

Also,  for  greater  flexibility,  let  us  use  a  parameter  q,  0  <  q  <  1,  in  the  seeond  faetor  so 
that 

^^{xln-x)  <ki}  =  qai  and  Pr{  >  ki}  =  (l-q)ai. 

This  provides  relief  from  the  practiee  of  splitting  the  tails  evenly  in  the  asymmetrie  ehi 
square  distribution. 

Eaeh  faetor  ean  be  pivoted  separately,  i.e., 

Pr{  X-Zgcr/ 4n  <  ju  <  X  +  <j I ^/n  }Pr{ ^SS / kj  <  <t  <  -.JSS / k^  }.  (6.17) 

The  eonfidenee  region  is  the  interseetion  of  the  two  events.  The  range  of  values  for  a  in 
the  seeond  faetor  is  used  in  the  interval  limits  of  the  first  faetor.  This  eapability  is 
programmed  into  the  funetion  NormCI  (x,  alph,  p,  q,  graph  =  T)  ,  a 
plotting  funetion  that  aeeepts  the  data  x  and  a,  where  1  -  a  is  the  eonfidenee  level.  If 
graph  =  F,  then  the  graph  is  not  drawn  and  the  return  is  the  set  of  vertiees  of  the 
trapezoid.  The  funetion  eontains  defaults  p  =  14,  q  =  14. 

6c.  Bias  Reduction 

The  estimators  for  p,  is  unbiased,  but  not  so  for  the  estimator  for  a.  In  this  latter 
ease,  we  ean  use  an  exaet  eorreetion  based  on  the  expeetation  the  maximum  likelihood 
estimator  for  a,  i.e.. 


E{s}  =  cr 


2  r(n/2) 

~n  r((n-l)/2) 


(6.18) 


and  the  exaet  bias  adjustment  is  obtained  from  (6.18). 

Note  1;  [Barndorff-Nielsen  and  Cox,  1994,  p.  187]  post  the  reeiproeal  of  the  third  factor 
in  (6.18). 
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Note  2:  If  one  is  using  the  unbiased  estimator  for  g^,  then  the  faetor  Vl/n  should  be 
replaeed  with  ^2/{n-l). 


6d,  Censored  Samples 

We  require  partial  derivatives  of  f  and  S,  (6.7)  for  use  in  (2.28),  (2.29),  and  (2.30). 

Inf  —  ju 

The  are  reeorded  here.  Begin  with  those  of  the  first  order,  using  z  = - —  and 


a 


w= 


ln(t)  -  ju 


a 


=  — T  ^)  =  (p(z)  =  z  f(x)/a 

xa  a  xa 

=  (z^-l)f(x)/G 
xa 

8^(0  =  —  (p(w) 
a 

Sa(t)  =  —  (p(w) 
a 

and  those  of  the  seeond  order, 

1)  =  (z^-1)  f(x)/a^ 
xa 

/^^(x)  =  ^y(p(z){l-z-3z^+z^}  =  2  f(x)[l -z-3z^  +  zVa^ 

xa 

=  z(z^  -  3)  f(x)/a^ 
xa 

w 

=  —(?iw) 
a 

1  7 

Saa  =  —W{W  -2)(p(w) 
a 

S^a  =  . 

a 

The  program  for  estimating  p,  and  a  is  I  test  ().  Its  use  is  illustrated  in  Seetion  3. 
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7.  Summary  and  Conclusions 

The  paper  offers  great  flexibility  to  those  practitioners  and  students  who  wish  to 
make  exploratory  computations  and  learn  the  behavior  of  models.  Those  possessing  a 
modicum  of  programming  capability  can  expand  the  choices  given  here  to  other  models 
and  situations  as  the  structure  offers  a  template  that  facilitates  additional  expansion. 
Much  remains  to  be  learned.  What  follows  is  a  listing  of  what  has  been  learned  thus  far. 
We  comment  on  the  items  as  presented  in  Section  3. 

Survival  distribution  selection  must  be  based  on  careful  modeling  and  good 
judgment.  This  area  of  statistical  work  has  the  disadvantage  that  the  collection  of  large 
quantities  of  data  requires  time,  often  unacceptable  amounts  of  time.  Yet  the  comparisons 
of  models  using  somewhat  small  data  sets  does  not  allow  one  to  examine  the  tail 
behavior.  The  sole  use  of  histograms  and  q-q  plots  does  not  suffice.  The  graphing  of  the 
other  support  functions,  especially  the  hazard  function  and  mean  residual  life,  can  aid  the 
practitioner  in  modeling  and  the  making  of  choices.  Such  is  illustrated  in  Sections  3a  and 
3b. 


The  graphing  of  the  behavior  of  planned  replacement  policies  can  be  quite  useful 
to  those  involved  in  the  generation  of  maintenance  policies.  The  problem  treated  in  3c 
deals  solely  with  the  relative  costs  of  planned  versus  random  replacements.  The  most 
interesting  point  is  the  flat  behavior  of  the  cost  curves  once  the  original  high  costs  fall 
away.  This  is  useful  to  the  maintenance  planner  because  the  replacement  model  is 
generically  a  very  simple  one,  one  that  does  not  consider  the  costs  of  arranging  to  make 
the  replacements  at  the  optimal  point  in  time.  Such  costs  are  likely  to  be  high  and  it  is 
comforting  to  know  that  the  optimal  point  is  contained  in  a  very  broad  window.  One  need 
not  give  ground  to  the  problems  of  performing  maintenance  at  inconvenient  times. 

Sections  3d,  3c,  and  3f  illustrate  the  implementation  of  the  general  materials 
presented  in  Section  2  for  the  Weibull,  Gamma,  and  Lognormal  distributions, 
respectively. 

There  are  four  areas  of  specialization  for  each; 

Computation  of  the  support  functions. 

Maximum  likelihood  parameter  estimation  for  complete  samples. 

Implementation  of  bias  reduction  techniques. 

Maximum  likelihood  estimation  for  right  censored  samples. 

Summary  of  each  is  made  in  turn. 

Most  software  systems  provide  the  capability  to  compute  density,  cumulative 
distribution,  and  quantile  functions.  The  extension  of  such  capabilities  to  the  set  of 
survival  distribution  support  functions  has  yet  to  occur.  The  three  cases  presented  here 
provide  illustration  as  to  what  must  be  done. 
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Maximum  likelihood  estimation  for  eomplete  samples  is  popular  and  well 
understood,  yet  they  are  hard  to  eompute.  The  algorithms  presented  here  have  high 
preeision.  The  statisties  texts  do  not  provide  mueh  in  the  way  of  support  for  joint 
eonhdenee  regions  for  multi-parameter  models.  The  author  believes  them  to  be  useful. 
They  bear  witness  to  Bellman’s  “eurse  of  dimensionality.”  It  is  quite  remarkable  how 
large  these  regions  are  for  samples  of  size  50  and  100. 

Maximum  likelihood  estimators  have  long  been  known  to  be  biased,  but  there 
appears  little  available  eommentary  on  the  importanee  of  this  bias.  The  issue  has  not 
reeeived  eonspieuous  attention.  The  method  implemented  here  was  taken  from  [Cox  and 
Hinkley];  it  is  a  first  order  approximation  based  on  the  marginal  distribution  of  the 
individual  parameters.  The  limited  eomputational  results  presented  are  not  very 
interesting;  but  the  stage  is  set  for  a  broadly  based  Monte  Carlo  simulation  study  geared 
to  learn  the  behavior  of  the  method.  At  present,  it  appears  that  it  may  be  useful  only  for 
rather  small  sample  sizes. 

Methods  for  managing  censored  data  may  be  available  in  the  Splida  and  Reliasoft 
systems.  But  these  are  opaque.  The  open  methods  that  appear  here  are  limited,  but  offer  a 
beginning.  The  functions  Itest  ( )  and  Golden  ( )  are  generic.  In  concept,  they  may  be 
applied  with  any  right-censored  data  set  modeled  with  any  distribution  that  is  smooth  in 
its  parameters,  and  any  finite  dimension  of  the  parameter  space.  However,  they  have  not 
been  broadly  tested.  In  fact,  their  use  was  not  competitive  with  the  bisection  search 
method  devised  for  the  Weibuh  distribution.  See  Section  4. 

The  idea  of  alternating  the  application  of  Newton-Raphson  iteration  with  a  golden 
section  search  may  be  new.  It  should  be  useful  whenever  the  log  likelihood  is  a  concave 
function  of  the  parameters.  It  appears  to  converge  in  a  reasonably  short  amount  of  time, 
in  spite  of  its  low  dependence  upon  the  quality  of  the  initializing  point. 

The  computation  of  derivatives  of  the  incomplete  gamma  function  came  up  in  the 
treatment  of  the  gamma  distribution  for  censored  samples.  A  method  of  high  precision 
was  developed  and  appears  in  Appendix  D.  It  is  suitable  for  parallel  processing.  Mostly 
the  values  are  precise  to  at  least  nine  decimal  places. 

The  computation  of  the  information  matrix  and  the  bias  reduction  term  is  quite 
difficult  for  censored  samples.  Indeed  the  calculation  depends  upon  the  type  of  censoring. 
Much  depends  on  the  information  matrix.  Also,  there  is  dependence  on  the  ratio  of  the 
number  of  uncensored  sample  to  the  total  sample  size.  The  quality  of  the  estimation 
results  will  deteriorate  as  the  number  of  uncensored  observations  becomes  small. 
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Locations  of  S-Plus  Listings  in  the  Appendices 


Appendix  A.  psi(),  psi.basQ,  gl(),  g2(),  g3(),  g4() 

Appendix  B.  weibulLmleQ,  weibull.mm(),  weibull.mle(),  weibull.estQ,  Newt.wei(), 
Init.wei(),  KapM(),  eonvtoweiQ 

Appendix  C.  Itest(),  GoldenQ,  init.gam(),  Llik.gam(),  Newt.gam(),  Lp.gamQ, 
Llik.logn(),  Newt.lognQ,  Lp.logn() 

Appendix  D.  Surv.gam(),  surv.gam(),  surva.gam(),  survaa.gam(),  transl(),  trans2(), 
trans3() 

Appendix  E.  gam.est(),  ellipseQ,  expl.lnomi(),  weil.gamQ,  wei2.gam(),  wei.cost(), 
NormCTO,  area.eomp(),  seg.eomp(),  sol.ptQ,  Newt.wei(),  exp.lnomi(), 
Ext.newtQ,  eost.eomp(),  Gsq.normO,  Gsq.gam(),  Gsq.wei() 
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Appendix  A 

Asymptotic  Expansions  for  the  Polygamma  Functions 

The  Psi  function  is  defined  as  the  derivative  of  the  log  gamma  function.  The 
recursion  formula  for  the  Gamma  function  translates  into  a  like  formula  for  the  Psi 
function.  Thus, 

\|/(z)  =  d(ln(r(z))/dz  and  \|/(z+l)  =  \|/(z)  +  1/z.  (A.l) 

Asymptotic  expansions  are  especially  advantageous  because  high  precision  is  available 
using  few  terms,  provided  z  is  large.  This  can  be  easily  exploited  in  computations 
involving  the  Polygamma  functions.  Thus,  if  z  +  r  is  sufficiently  large  to  obtain 


appropriate  accuracy,  then  use 

v|/(z)  =  \|/(z  +  r)  -  1/z  -  l/(z  +  1)  -  ...  -  l/(z  +  r  -  1)  (A.2) 

v|/'(z)  =  \|/'(z  +  r)  +  1/z^  +  l/(z  +  1)^  +  . . .  +  l/(z  +  r  -  1)^  (A.3) 

and  so  on.  It  remains  to  record  the  expansions  for  large  z,  [Abramowitz  and  Stegun] 

v|/(z)  =  ln(z)  -  l/2z  -  1/  12z^  +  l/120z^  -  l/252z^  +  . . .  (A.4) 

V'(z)  =  1/z  +  l/2z^  +  l/6z^  -  l/30z^  +  l/42z’  -  l/30z‘^  +  . . .  (A.5) 

\|f"(z)  =  -l[l/z^  +  1/z^  +  l/2z^  -  l/6z^  +  l/6z^  -  3/lOz^'’  +] 

V'''(z)  =  Hz  +  Hz  +  Hz  -  1/z’  +  4/3z‘^  -  3/z“  + 


V^"\z)  =  (-1)"+'  [(n-l)!/z"  +  n!/2z“+^  +  EB2k  (2k  +  n  -  l)!/(2k)!z’'^+"](A.6) 
and  {B2k}  are  the  Bernoulli  numbers.  See  [Abramowicz  and  Stegun]. 

The  error  in  these  expansions  is  sized  by  the  first  term  ignored.  They  work  best 
for  z  large,  say  z  >  10.  For  smaller  values  of  z  one  should  make  adjustments  based  on  the 
recursive  formula  for  the  gamma  function. 

r(z  +  a)  =  (z  +  a  - 1)  •  •  •  (z  +  l)(z)r(z) 

'k(z  +  a)=  +  'P(z)  (A.7) 

0  Z  +  j 

/  ^  if 

'P'(z  +  a)  =  -Z  -  +'P'(z) 

0  Iz  +  jj 
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\|f"(z+a) 


2Z( 


z  +  J 


V"'(z+a)  =  -6|;( 

0 


Z  +  J 


For  these  reasons,  two  programs  have  been  written:  psi.bas  (z)  takes  a  veetor 
argument  z  and  returns  three  rows  of  output,  the  funetion  and  its  first  two  derivatives.  For 
small  values  of  z,  one  ean  ehoose  an  integer  a  and  call  psi  (z,a)  to  achieve  precision  based 
on  (A.7).  The  default  value  of  a  is  10.  Adjustments  are  made  for  all  components  of  z  less 
than  a. 


On  occasion  we  require  the  derivatives  of  the  gamma  function.  These  are 
obtainable  from  the  psi  function  and  its  derivatives.  The  first  four  are  as  follows. 


r(z)  =  r(z)v(z) 

r"(z)  =  r(z){V(z)  +  V(z)} 

r"'(z)  =  r(z){v"(z)  +  3v'(zMz)  +  V(z)}  (A.8) 

r"'(z)  =  r(z){\|/"’(z)  +  4\|/"(z)\|/(z)  +  6\|/'(z)V(z)  +  3[\|/'(z)]^  + V(z)}. 


Useful  relationships 

W{x)e-^dx  =  and  ^ 


7=1 


rO)(v)vpi'-./i(v) 


(r-J), 


f(«)=>f(i)+^;:T;  'f"(n)='i"(i)-z;;|^; 


S-Plus  Codes 

psi 

function(z,  a  =  10){ 

#  finame  is  psi;  programmer;  R.  Read 

#  derivative  log  gamma  and  two  subsequent  derivatives 

#  Uses  recursive  formula  for  min(z)  <  a  and  the 

#  psi.bas  function  for  large  z  (asymptotic  expansion) 

a  <-  floor(a) 
s  <-  a  -  floor(min(z)) 
j  <-  l:s 
if(s  >  0)  { 

tempo  <-  apply(l/outer(z  -  1,  j,  "+"),  1,  sum) 
tempi  <-  apply(l/outer(z  -  1,  j,  "+")^2,  1,  sum) 
temp2  <-  apply(l/outer(z  -  1,  j,  "+")^3,  1,  sum)  } 

else  { 

tempo  <-  tempi  <-  temp2  <-  0} 
if(s  <=  0) 

s  <-  0 
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out  <-  psi.bas(z  +  s) 
out[l,  ]<-out[l,  ]- tempo 
out[2,  ]  <- out[2,  ]  + tempi 
out[3,  ]  <-  out[3,  ]  -  2  *  temp2 
out} 


psi.bas 

funotion(z){ 

#  fname  is  psi.bas;  programmer:  R.  Read 

#  asymptotie  expansion  for  first  three  derivatives 

#  of  the  log  gamma  funetion;  z  may  be  a  veetor. 

K  <-  length(z) 

eoefO  <-  c(-2, -12,  120,  -252) 

mO  <-  matrix(rep(eoefO,  rep(K,  4)),  neol  =  4) 

zO  <-  e(z,  z^2,  zM,  z^6) 

eoefl  <-  e(l,2,  6,  -30,  42,  -30) 

ml  <-  matrix(rep(eoef  1 ,  rep(K,  6)),  neol  =  6) 

zl  <-  o(z,  z^2,  z^3,  z^5,  z^7,  z^9) 

eoefZ  <-  o(-l,  -1,  -2,  6,  -6) 

m2  <-  matrix(rep(ooef2,  rep(K,  5)),  neol  =  5) 

z2  <-  o(z^2,  z^3,  z^4,  z^6,  z^8) 

out  <-  log(z)  +  apply} l/(m0  *  zO),  1,  sum) 

out  <-  rbind(out,  apply(l/(ml  *  zl),  1,  sum)) 

out  <-  rbind(out,  apply(l/(m2  *  z2),  1,  sum)) 

out} 


gl 

funotion(z){ 

#  fname  is  gl 

#  first  derivative  of  the  gamma  funetion 

out  <-  gamma(z)  *  psi(z)[l,  ] 
return}  out)} 


g2 

funotion}z){ 

#  fname  is  g2 

#  seeond  derivative  of  the  gamma  funetion 

tmp  <-  psi}z) 

out  <-  gamma}z)  *  }tmp[2,  ]  +  tmp[l,  ]^2) 
return}out)} 


g3 

funotion}z){ 

#  fname  is  g3 

#  third  derivative  of  the  gamma  funetion 

tmp  <-  psi}z) 
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out  <-  gamma(z)  *  (tmp[3,  ]  +  3  *  tmp[2,  ]  *  tmp[l,  ]  +  tmp[l,  ]^3) 
return(out)} 


g4 

function(z){ 

#  fname  is  g4 

#  fourth  derivative  of  the  gamma  function 

tmp  <-  psi(z) 

out  <-  gamma(z)  *  (tmp[4,  ]  +  4  *  tmp[3,  ]  *  tmp[l,  ] 

+  6  *  tmp[2,  ]  *  tmp[l,  ]^2  +  4  *  tmp[2,  ]^2  +  tmp[l,  ]M) 
return(out)} 
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Appendix  B 


Analysis  for  Computational  Support  of  the  Weibull  Distribution 

This  appendix  contains  the  numerical  analysis  and  S-Plus  code  for  executing  the 
support  functions  developed  in  Chapter  4.  Fundamental  support  functions  for  the  Weibull 
distribution  are  primitives  in  the  S-Plus  system.  Capabilities  for  computing  the  psi  and 
gamma  functions  and  their  derivatives  are  contained  in  Appendix  A.  The  fact  that 
Y  =  (X/p)“  has  an  Exp(l)  distribution  when  X  ~  Weibull(a,  P)  has  been  exploited  and 
the  properties  of  the  Exponential  distributions  are  familiar. 

It  remains  to  describe  the  computational  techniques  used  to  compute  maximum 
likelihood  estimates  for  the  shape  and  scale  parameters  (a,  P)  under  complete  and 
censored  sampling. 

Complete  Samples 

The  maximum  likelihood  estimates  are  computed  using  the  natural  iteration 
function  described  in  Chapter  4.  The  method  of  moments  estimates  are  used  for 
initialization.  The  Newton-Raphson  Iteration  (B.l)  is  used  to  find  these  latter  estimates. 
The  functions  utilized  are 

weibull.mleO;  weibull.mmQ;  weibull. est(); 

the  last  of  these  functions  calls  the  previous  two  and  displays  both  kinds  of  estimates. 

Their  use  is  illustrated  in  Section  3  and  the  S-plus  listings  appear  in  the  S-Plus  section  of 
this  appendix. 

Right-Censored  Samples 

The  formulae  in  Section  4d  require  computational  development.  The  transform 
the  extreme  value  distribution  enables  one  to  execute  a  one-dimensional  search  for  the 
parameter  b,  see  (4.24)  and  (4.28).  However,  it  appears  that  the  Newton-Raphson  scheme 
undergoes  considerable  oscillation  when  applied  in  this  setting  and  some  modifications 
are  in  order.  The  basic  idea  is  to  use  Newton-Raphson  until  the  function  h  of  (4.28) 
oscillates  and  then  switch  to  a  bisection  search. 

A  little  more  detail  can  be  useful.  The  basic  iteration  is  to  use 

bn+l  =  b„-h(bn)/h'(b„)  (B.l) 

and  record  the  two  most  recent  pairs,  call  them  (bi,  hi)  and  (b2,  h2).  At  each  step  check 
the  sign  of  the  product  hixh2;  when  it  turns  negative,  change  to  the  bisection  search.  I.e., 


50 


set  b  =  (bi  +  b2)/2;  compute  h(b); 

if  hxhi  <  0  set  h2  =  handb2  =  b;  otherwise  set  hi  =  handbi  =  b;  (B.2) 

repeat  until  no  change  in  b,  or  h  «0,  or  both. 

The  S-Plus  function  that  does  this  is  called  Newt.wei(). 

The  initialization  exploits  Equation  (4.30).  There  we  have  a  straight-line 
relationship  between  v  and  y,  which  can  be  expressed  using  a  least  squares  fit.  Then  set 
b  =  1/slope  and  u  =  x-intercept.  The  conversion  of  {vi}  to  {yi}  is  accomplished  using 
the  Kaplan-Meier  estimator  for  the  survivor  function.  It  was  decided  to  write  a  simple 
code  to  accomplish  this. 

The  basic  Kaplan-Meier  estimation  has  the  following  rules.  Pool  and  order  the 
uncensored  and  censored  times,  call  them  {x*}  .  Let  di  be  the  number  that  died  at  x*,  i.e., 
exclusive  of  those  that  were  censored  at  that  value.  Set  no  =  n  and  ni  =  ni.i  -  di  for  all  of 
the  unique  values  of  {x*} .  These  are  the  number  at  risk  values.  Then  the  estimated 
survivor  function  is,  at  the  data  points. 


s(x;) 


-d^ 


(B.3) 


These  details  are  executed  by  the  functions  init.wei()  and  KapM();  listings  below. 

S-Plus  Codes 


weibull.mle 

function(x,  aO  =  1,  ep  =  0.0001){ 

#  finame  is  weibull.mle 

#  natural  iteration  function;  aO  is  initial  shape  parameter 

#  output  is  the  pair  (shape,  scale) 

Ixb  <-  mean(log(x)) 
a  <-  aO 
repeat  { 

ainv  <-  mean(log(x)  *  x^a)/mean(x^a)  -  Ixb 
a  <-  1/ainv 
if(abs(a  -  aO)  <  ep) 
break 

aO  <-  a 

} 

b  <-  (mean(x^a))^(l/a) 

out  <-  c(a,  b) 

out} 

weibull.mm 

function(x,  aO  =  1,  ep  =  0.0001){ 
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#  fname  is  weibull.mm 

#  returns  the  method  of  moments  estimator  for  the  shape  (a)  and  seale  (b) 

#  parameters  of  the  weibull  distribution.  The  initialization  is  aO.  The  data  are  x. 

xb  <-  mean(x) 
s2  <-  var(x) 
a  <-  aO 
repeat  { 

g  <-  lgamma(l  +  2/a)  -  2  *  lgamma(l  +  1/a)  -  log(l  +  s2/xb^2) 
gp  <-  (psi(l  +  2/a)[l]  -  psi(l  +  l/a)[l])  *  (-2/a^2) 
a  <-  aO  -  g/gp 
if(a  <  0) 

a  <-0.1 

if(abs(a  -  aO)  <  ep) 
break 
aO  <-  a  } 

b  <-  xb/gamma(l  +  1/a) 

out  <-  o(a,  b) 

out} 

weibull.mle 

funotion(x,  aO  =  1,  ep  =  0.0001) 

{#  fname  is  weibull.mle 

#  natural  iteration  funetion;  aO  is  initial  shape  parameter 

#  output  is  the  pair  (shape,  scale) 

Ixb  <-  mean(log(x)) 
a  <-  aO 
repeat  { 

ainv  <-  mean(log(x)  *  x^a)/mean(x^a)  -  Ixb 
a  <-  1/ainv 
if(abs(a  -  aO)  <  ep) 
break 
aO  <-  a } 

b  <-  (mean(x^a))^(l/a) 

out  <-  c(a,  b) 

out} 

weibull.  est 
function(x){ 

#  fname  is  weibull. est 

#  Input  is  a  random  sample  from  the  Weibull  distribution. 

#  Output  has  five  components:  shape  and  scale  parameter  estimates 

#  using  max  likelihood,  sample  size,  shape  and  scale  parameter 

#  estimates  using  method  of  moments 

n  <-  length(x) 
mm  <-  weibull. mm(x) 
ml  <-  weibull. mle(x) 
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out  <-  c(ml,  n,  mm) 

names(out)  <-  c("shape:mle",  "scaleimle",  "samp  size",  "shapeimm",  "scale:mm") 
out} 

Newt.wei 

function(x,  t,  param,  ep  =  0.000 1){ 

#  fname  is  Newt.wei 

#  Newton-Raphson  method,  then  biseetion  search  used  to  find  mle 

#  for  the  Weibull  distribution  using  the  extreme  value  distribution  technique. 

#  Data  input  is  (x,  t);  param  initialization  is  (u,  b),  b>0 

#  output  is  (u,  b,  alph,  beta),  r  is  cardinality  of  uncensored  set 

r  <-  length(x) 

V  <-  log(x) 
vb  <-  mean(v) 

It  <-  log(t) 
b  <-  param[2] 
u  <-  param[l] 
vv  <-  c(v.  It) 
b2  <-bl  <-b 
hi  <-  h2  <-  0 

j<-l 

flag  <-  1 
repeat  { 

if(flag  ==  1)  { 

D  <-  sum(exp(vv/b)) 

N  <-  sum(vv^2  *  exp(vv/b)) 

A  <-  sum(vv  *  exp(vv/b)) 
h  <-  A/D  -  b  -  vb 

bp  <-  -  N/(D  *  b^2)  +  A^2/(D  *  b)^2  -  1 
b<-bl  -h/hp 
if(b  <=  0) 

b<-0.01 

if(max(abs(h),  abs(b  -  bl))  <  ep) 
break 

if(j/2  !=  round(j/2))  { 
bl  <-b 
hi  <-h} 

if(j/2  ==  round(j/2))  { 
b2  <-b 
h2<-h} 

j  <-j  +  1 

if(sign(h2  *  hi)  <  0) 
flag  <-  0} 

ifG==25) 

break 
if(flag  ==  0)  { 
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b  <-  (bl  +  b2)/2 
D  <-  sum(exp(vv/b)) 

A  <-  sum(vv  *  exp(vv/b)) 
h  <-  A/D  -  b  -  vb 
if(sign(h  *  hi  <  0))  { 
h2  <-h 
b2  <-b} 

if(sign(h  *  h2  <  0))  { 
hi  <-h 
bl  <-b} 

if(max(abs(h),  abs(b  -  bl))  <  ep) 
break 

j  <-j  +  1 

ifO  ==  50) 

break}  } 

u  <-  b  *  log(D/r) 
alph  <-  l^ 
beta  <-  exp(u) 

out  <-  e(u,  b,  alph,  beta,  flag) 

names(out)  <-  e("u","b",”alph","beta","flag") 

return(out)} 


init.wei 
funetion(x,  t){ 

#  fname  is  init.wei 

#  provides  initial  estimates  for  the  extreme  value  distribution 

#  i.e.,  y  =  log(weibull),  parameters  (u,  b).  The  pair  (x,  t) 

#  represents  n  observations  (duplieate  values  required)  where 

#  the 't'  values  are  the  eensoring  values.  Output  is  a  four  veetor; 

#  first  two  are  (u,  b)  and  the  last  two  are  (alpha,  beta). 

#  Isfit  to  the  log  survival  fne  teehnique  is  utilized. 

S  <-  KapM(x,  t) 
w  <-  log(  -  log(S)) 
y  <-  sort(log(e(x,  t))) 
yy  <-  unique(y) 

yb  <-  mean(yy)  #  w  <-  lsfit(yy,  S)$eoef 

slope  <-  sum(w  *  (yy  -  yb))/sum((yy  -  yb)^2) 

intere  <-  mean(w)  -  slope  *  yb 

b  <-  1/slope 

u  <-  -  b  *  intere 

alph  <-  slope 

beta  <-  exp(u) 

names(out)  <-  e("u",  "b",  "alph",  "beta") 
out  <-  o(u,  b,  alph,  beta) 
return(out)} 

KapM 
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function(x,  t){ 

#  fname  is  KapM 

#  produces  the  Kaplan-Meier  estimate  of  the  survivor  funetion  for 

#  data  o(x,  t)  where  x  are  the  aetual  death  times  and  t  are  the 

#  eensored  values. 

y  <-  sort(e(x,  t)) 
yy  <-  unique(y) 

dd  <-  d  <-  tab  <-  table(y)  #  number  of  deaths  at  y 
n  <-  length(y) 
k  <-  length(yy) 

nr  <-  n  -  tab  #  initial  number  at  risk 

tt  <-  unique(t) 

ind  <-  (1  ;k)[tt  ==  yy] 

tabt  <-  table(tt) 

kk  <-  length(ind) 

if(kk  >  0)  { 

for(j  in  l:kk) 

dd[indD]]  <-  d[indD]]  -  tabtQ]} 

SS  <-  (nr  -  dd)/nr 
S  <-  eumprod(SS) 
return(S)} 
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Appendix  C 


Implementation  of  the  General  Censored  Data  Estimation  Scheme 

The  iterative  estimation  scheme  described  in  Section  2  has  been  implemented  for 
two  out  of  three  distributions,  namely  the  gamma  and  lognormal  distributions.  The  case 
of  the  Weibull  distribution  can  be  managed  more  directly  as  described  in  Section  4  and 
Appendix  B. 

Much  of  the  structure  can  be  treated  generally.  Because  S-Plus  functions  have  the 
capability  to  accept  other  functions  as  input,  it  is  possible  to  write  generic  code  for  the 
process.  Such  is  exploited  for  the  values  of  the  log  likelihood,  the  golden  section  search, 
and  the  change  of  gradient  direction  effected  by  a  single  iteration  of  the  Newton-Raphson 
scheme.  The  name  Itest  is  short  for  iterative  estimation. 

The  basic  call  is 

Paramest  <- Itest(  X,  t,  thO,  Newt,  Llik,  ep  =  .0001), 

where  x  is  the  set  of  uncensored  survival  times;  t  is  the  set  of  right-censored  times 
realized;  thO  is  the  initialization  point  in  the  parameter  space;  Newt  is  a  generic  function 
that  produces  an  updated  value  for  the  parameter;  and  Llik  is  a  generic  function  that 
computes  the  log  likelihood  for  the  targeted  distribution  family.  The  value  of  ep  is  used  to 
control  the  precision  of  the  estimate.  Of  course,  the  inputs  to  Newt  and  Llik  must  be 
generated  by  the  parent  function. 

The  function  also  calls  a  golden  section  search  function.  Its  role  is  to  seek  the  best 
value  for  th  along  the  line  segment  connecting  th  previous  and  the  output  of  Newt. 
Having  made  that  selection,  the  parent  program  calls  Newt  to  find  a  new  direction  for 
search.  The  output  has  the  structure 

Paramest  =  th[l],  th[2],  N.iter, 

where  N.iter  is  the  number  of  cycles  through  Newt. 

S-Plus  Codes 

Itest 

function(x,  t,  thO,  Newt,  Llik,  ep  =  0.0001,  N.iter  =  50) 

{#  fname  is  Itest 

#  alternates  the  use  of  golden  section  and  Newton-Raphson 

#  to  find  2-D  mle's.  The  input  N.iter  puts  a  cap  on  the 

#  number  of  iterations  of  Newton-Raphson. 

recth  <-  NULL 
k<-  1 

repeat  {if(k  ==  50)  break 
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thl  <-  Newt(x,  t,  thO) 
if(max(abs(thO  -  thl))  <  ep) 
break 

recth  <-  rbind(recth,  c(thO,  thl)) 
k  <-  k  +  1 

thO  <-  Golden(x,  t,  thO,  thl,  Llik)[[l]]} 
out  <-  c(recth[k  -  1,  3:4],  round(k)) 
if(k  ==  N.iter) 

out  <-  c("No  convergence") 
return(out)} 


Golden 

function(x,  t,  thO,  thl,  Llik,  epl  =  0.5) 

{#  fname  is  Golden 

#  performs  the  golden  section  search  along  the  direction 

#  thO  to  thl .  The  function  returns  the  value  of  th  that 

#  maximizes  Llik  in  that  direction.  The  return  is  the 

#  new  th  and  its  (proportional)  distance  from  thO. 

#  This  function  is  called  by  Itest,  from  which  it  gets  ep 

pi  <-  c(0.618,  0.382) 

p2  <-  rev(pl)  #thmit  <-  cbmd(th0,  thl) 

to  <-  thO 

lenO  <-  sqrt(sum((thl  -  th0)^2)) 

j  <-  1 

repeat  {j  <-j  +  1 

th4  <-  pi  *  thO  +  p2  *  thl 
th6  <-  p2  *  thO  +  pi  *  thl 
th  <-  cbmd(th0,  th4,  th6,  thl) 

L  <-  Lhk(x,  t,  th) 

ifG>5){ 

ep2  <-  (max(L)  -  mm(L))/20 
ep2  <-  max(ep2,  ep) 
epl  <-  mm(epl,  ep2)} 
ifO  ==  100) 

break 

rnk  <-  rank(L)  #  print(L) 
if(mk[3]  ==  4) 
thO  <-  th4 
if(mk[2]  ==  4) 

thl  <- th6 
if(mk[4]  ==  4) 
thO  <-  th4 
if(mk[l]  ==  4) 

thl  <- th6 

if(max(abs(thO  -  thl))  <  epl)  break  #  prmt(c(th0,  thl))} 
lenl  <-  sqrt(sum((thl  - 10)^2)) 
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out  <-  list(th  =  thl,  prop  =  lenl/lenO) 
ifO  ==  100)  { 

out  <-  c("Too  many  Golden  iterations")} 
return(out)} 

Input  functions  when  using  the  Gamma  distribution. 

Begin  with  the  computation  of  the  initializing  point. 

init.gam 

function(x) 

{#  fname  is  init.gam 

#  initializes  the  estimates  using  the  method  of  moments  applied  to 

#  the  uncensored  part  of  the  data 
#x<-  dat[dat[,  2]==  1,  1] 

xb  <-  mean(x) 
s2  <-  var(x) 
lam  <-  xb/s2 
alph  <-  lam  *  xb 
out  <-  c(alph,  lam) 
out} 

Llik.gam 
function(x,  t,  th) 

{#  fname  is  Llik.gam 

#  computes  the  log  likelihood  for  censored  gamma  data 

#  X  in  the  uncensored  life  lengths,  t  is  the  censored  set 

#  th  is  a  matrix  of  parameter  values; 

#  first  row  is  alph,  second  is  lam 

if(length(th)  ==  2)  th  <-  matrix(th,  2,  1) 
k  <-  ncol(th) 
r  <-  length(x) 

Ig  <-  lgamma(th[l,  ]) 
lx  <-  sum(log(x)) 
sx  <-  sum(x) 

L  <-  rep(0,  k) 
for(i  in  1  :k)  { 

S  <-  1  -  pgamma(th[2,  i]  *  t,  th[l,  i],  1) 

L[i]  <-  r  *  (th[l,  i]  *  log(th[2,  i])  -  r  *  lg[i])  + 

(th[l,  i]  *  -1)  *  lx  -  th[2,  i]  *  sx  +  sum(log(S))} 

return(L)} 

Newt.gam 
function(x,  t,  th) 

{#  fname  is  Newt.gam 

#  program  executes  a  single  iterative  update  of  th  (theta) 
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#  for  the  mle  estimation  of  th  (alph,  lam)  for  eensored  gamma 

#  distribution  data;  x  is  the  set  of  observed  life  times  and  t 

#  the  set  of  eensored  (right)  times. 

alph  <-  th[l] 
lam  <-  th[2] 
ps  <-  psi(alph) 
r  <-  length(x) 

S  <-  Surv.gam(alph,  lam  *  t,  20,  ep  =  le-005) 

La  <-  r  *  (log(lam)  -  ps[l])  +  sum(log(x))  +  sum(S[2,  ]/S[l,  ]) 
f  <-  dgamma(lam  *  t,  alph,  1) 

Slam  <-  - 1  *  f 

LI  <-  (r  *  alph)/lam  -  sum(x)  +  sum(Slam/S[l,  ]) 

Laa<-  -  r  *  ps[2]  +  sum(S[3,  ]/S[l,  ]  -  (S[2,  ]/S[l,  ])^2) 
Salphlam  <-  - 1  *  f  *  (log(lam  *  t)  -  ps[l]) 

Slamlam  <-  - 1^2  *  f  *  ((alph  -  l)/(lam  *  t)  -  1) 

Lai  <-  r/lam  +  sum(Salphlam/S[l,  ]  -  (S[2,  ]*  Slam)/S[l]^2) 

Lll  <-  ( -  r  *  alph)/lam^2  +  sum(Slamlam/S[l,  ]  -  (Slam/S[l,  ])^2) 
Lt  <-  c(La,  LI) 

H  <-  matrix(c(Laa,  Lai,  Lai,  Lll),  2,  2) 
thO  <-  th 

th  <-  thO  -  solve(H,  Lt)  #  print(th) 
if(th[l]  <  0) 

th[l]<-0.01 
if(th[2]  <  0) 

th[2]  <-  0.01 

dist  <-  max(abs(th  -  thO))  #  print(dist) 
retum(th)} 

The  funetion  Surv .  gam  ( )  is  deseribed  in  Appendix  D. 

Lp.gam 

funetion(x,  t,  thO) 

{#  fname  is  Lp.gam 

#  ereates  the  gradient  veetor  Lp  and  the  Hessian  matrix  H 

#  for  the  eensored  gamma  family  sampling 

alph  <-  th[l] 
lam  <-  th[2] 

S  <-  Surv.gam(alph,  lam  *  t,  aO  =  20) 
xml  <-  (log(lam  *  x)  -  psi(alph)[l,  ])/lam 
xni2  <-  alph/lam  -  x 
tml<-(S[2,  ]-psi(alph)[l,  ])/S[l,  ] 
tm2  <-  ( - 1  *  dgamma(lam  *  t,  alph))/S[l,  ] 

Lpl  <-  sum(xml)  +  sum(tml) 

Lp2  <-  sum(xm2)  +  sum(tni2) 

HI  1  <-  sum(xml^2)  +  sum(tml^2) 

H22  <-  sum(xm2^2)  +  sum(tm2^2) 
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H12  <-  sum(xml  *  xm2)  +  sum(tml  *  tni2) 

out  <-  list(Lp  =  c(Lpl,  Lp2),  H  =  matrix(c(Hl  1,  H12,  H12,  H22),  2,  2)) 
out} 

Input  functions  when  using  the  Lognormal  distribution. 

Initialize  with  mean(log(x))  and  stdev(log(x)) 

Llik.logn 
function(x,  t,  th) 

{#  fname  is  Llik.logn 

#  Computes  the  log  likelihood  for  censored  log  normal  data 

#  X  is  the  set  of  uncensored  life  lengths,  t  is  the  eensored  set. 

#  th  is  the  matrix  of  parameter  values;  first  row  is  mu, 

#  second  is  sigma.  The  output  is  a  set  of  log  likelihood  values. 

if(length(th)  ==  2)  th  <-  matrix(th,  2,  1) 
k  <-  ncol(th) 
r  <-  length(x) 
d  <-  length} t) 

zl  <- outer(log(x),  th[l,  ],  "-")/th[2,  ] 
z2  <-  outer(x,  th[2,  ],  "*") 
ec  <-  log(2  *  pi) 

LI  <-  apply}  -  log}z2)  -  0.5  *  cc  -  0.5  *  zl,  2,  sum) 
w  <- outer}log}t),  th[l,  ],  "-")/th[2,  ] 

Sn  <-  matrix}!  -  pnorm}w),  d,  k) 

L2  <-  apply}log}Sn),  2,  sum) 

L  <-  LI  +  L2 
return}L)} 

Newt.logn 
funetion}x,  t,  th) 

{#  fname  is  Newt.logn 

#  program  executes  a  single  iterative  update  of  th  }theta) 

#  for  the  mle  estimation  of  th  }mu,  sig)  for  eensored  lognormal 

#  distribution  data;  x  is  the  set  of  observed  life  times  and  t 

#  the  set  of  }right)  eensored  times. 

mu  <-  th[l] 
sig  <-  th[2] 

out  <-  Lp.logn}x,  t,  th) 

Lt  <-  out[[l]] 

H  <-  out[[2]] 
thO  <-  th 

th  <-  thO  -  solve}H,  Lt)  #  print}th) 
if}th[2]  <  0) 

th[2]  <-  0.01 

dist  <-  max}abs}th  -  thO))  #  print} dist) 
retum}th)} 
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Lp.logn 

function(x,  t,  thO) 

{#  fname  is  Lp.logn 

#  Creates  the  gradient  vector  Lp  and  the  Hessian  matrix  H 

#  for  the  censored  lognormal  family  sampling. 

#  The  output  is  a  list;  Lp  and  H 

mu  <-  thO[l] 
sig  <-  th0[2] 
z  <-  (log(x)  -  mu)/sig 
w  <-  (log(t)  -  mu)/sig 
fw  <-  dnorm(w) 

Sw  <-  1  -  pnorm(w) 

Lpl  <-  sum(z)  +  sum(fw/(sig  *  Sw)) 

Lp2  <-  sum((z^2  -  l)/sig)  +  sum((w  *  fw)/(sig  *  Sw)) 

HI  1  <-  sum((z^2  -  l)/sig^2  -  z^2)  +  sum((w  *  fw)/(sig^2  *  Sw)) 

H11<-H11-  sum(fw/(sig  *  Sw)^2) 

H12  <-  sum((z  *  (z^2  -  3))/sig^2  -  (z  *  (z^2  -  l))/sig) 

H12<-H12  -  sum(((w^2  -  1)  *  fw)/(sig^2  *  Sw)) 

H22  <-  sum(2  *  (1  -  z  -  3  *  z^2  +  z^3)  -  ((z^2  -  l)/sig)^2) 

H22  <-  H22  +  sum((w  *  (w^2  -  2)  *  fw)/(sig^2  *  Sw)  -  ((w  *  fw)/(sig  *  Sw))^2) 
out  <-  list(Lp  =  c(Lpl,  Lp2),  H  =  matrix(c(Hl  1,  H12,  H12,  H22),  2,  2)) 
return(out)} 
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Appendix  D 


Derivatives  of  the  Incomplete  Gamma  Function 

The  function  Surv.gamO  is  needed  in  Section  5  for  use  in  the  parameter 
estimation  programs  for  censored  gamma  data.  The  analysis  is  too  lengthy  to  be  included 
there  and  the  development  may  have  general  interest.  Basically  we  need  to  compute  the 
first  and  second  derivatives  of  the  Incomplete  Gamma  function  with  respect  to  the  shape 
parameter. 

The  formulas  used  are  taken  from  [Abramowicz  and  Stegun]  and  their  notation  is 
adopted.  We  focus  on  [Abramowicz  and  Stegun,  (6.5.1),  6.5.4),  and  (6.5.29)] 

P(a,  x)  =  — — f  y^^  'e^^dy  (D-1) 

r(a)Jo 


r-(a,x)  =  x-P(a,x)=  (D.2) 

The  power  series  in  (D.2)  plays  a  key  role.  It  will  be  used  for  values  of  x  smaller  than  a. 
With  this  understanding,  let  us  develop  a  crude,  but  useful,  bound  for  the  remainder  after 
n  terms. 


S" 


x^ 


''r(a  +  j  +  l)  r(a  +  n) 


1  CO 

-l-nt  ■^j=n 


X 


^  (a  “1“  n)  03  x^ 

(a  +  n)---(a  +  j)  r(a  +  n) (a  +  n)^ 


X 


a  +  n 


< 


r(a  +  n)'^'^"°  a  +  n  r(a  +  n)a  +  n-x  r(a  +  n) 


(D.3) 


This  inequality  allows  one  to  choose  n  given  a,  x  with  x  <  a  so  that  the  power  series  can 
be  made  as  precise  as  desired. 

The  issue  of  x  larger  than  a  is  managed  by  the  recursive  formula  [Abramowicz 
and  Stegun,  (6.5.21)] 


P(a+l,x)  =  P(a,  x)  - -, 

r(a  +  l) 

which,  upon  repeated  use,  yields 


(D.4) 


Inc(r)  =  P(a-r,x)  -  P(a,x)  =  e  ""  V  - . 

r(a  +  l-7) 


(D.5) 
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Viewing  the  goal  as  the  eomputation  of  P(a-r,  x),  one  may  always  ehoose  r  so  that  x  <  a 
and  then  the  power  series  portion  ean  finish  the  job.  Sueh  is  the  taetieal  plan.  Both  (D.2) 
and  (D.5)  may  be  differentiated  termwise. 


(D.6) 


r(a  +  7  +  l) 


['P^(a  +  7  +  l)-'P'(«  +  7  +  l)] 


(D.7) 


Inc^ir)  =  e  "X  ■=iV7~r — + 1  “ 7)] 


r(a  +  l-7) 


(D.8) 


a- 1 

\r 


Inc^Sr)  =  e-^Y<  ,  .  I[ln(x)-fi^(a  +  l-7)f +  •  (D.9) 

^-'r(a  +  l-7) 


Using  (D.l)  and  (D.2)  we  ean  express  the  survival  funetion  and  derivatives  of  the 
standard  gamma  random  variable. 


S(a,  x)  =  1-x"y>,  x)  =  1-e^^XIo  fU - ^ 

r(a  +  7  +  1) 


a+ j 

Sa(a,x)  =  -e-^^J^^--^^_[ln(x)-'P(a  +  l  +  7)]  (D.ll) 


Saa(a,x)  ,^_^{[ln(x)-'P(a  +  l  +  7)f-'P'(«  +  l  +  7)}  (D.12) 

^  r(a  +  l  +  j) 

Our  goal  is  to  orehestrate  these  formulas  into  a  eomputational  paekage  usable  for  the 
censored  sampling  likelihood  equations.  Moreover,  the  package  should  accept  vector 
input  for  x.  The  smaller  values  of  x  are  easily  managed  using  the  power  series.  For  other 
values  it  is  wise  to  use  the  tactic  of  shifting  to  larger  values  of  the  parameter  a  using 
(D.5),  (D.8),  and  (D.9).  Accordingly,  we  partition  the  (a,  x)  plane  into  three  horizontal 
strips 


vl;  0  <  X  <  8;  v2:  8  <  x  <  ao;  v3:  ao  <  x, 

and  aO  can  be  selected  for  power  series  truncation  error  control.  Our  programs  use 
aO  =  20.  Further,  each  strip  is  partitioned  selectively  for  purposes  of  invoking  the  tactic 
of  shifting  to  larger  values  of  the  parameter  a.  The  rule  in  place  is 

vl;  a  <  ao;  v2;  a  <  2ao;  v3:  a  <  3ao. 
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The  main  program  that  executes  these  rules  is  Surv.gam().  Sixteen  terms  are  used 
in  the  truncate  power  series.  For  these,  Equations  (D.IO),  (D.ll),  and  (D.12)  are 
computed  using  the  functions 

surv.gam():  surva.gam();  survaa.gamQ,  respectively. 

When  it  is  appropriate  to  shift  to  larger  values  of  a,  we  invoke  (D.5),  (D.8),  and 
(D.9).  These  functions  are  called 

translQ;  trans2();  trans3(),  respectively. 

The  output  of  Surv.gam  is  a  three-row  matrix  and  number  of  columns  equal  to  the 
cardinality  of  the  data  x.  The  three  rows  contain  the  zero*^,  first  and  second  derivatives. 

S-Plus  Codes 

Surv.gam 

function(a,  x,  aO  =  20,  ep  =  le-009) 

{#  fname  is  Surv.gam 

#  Integrates  the  programs  surva.gam,  survaa.gam 

#  and  trans2,  trans3  to  comput  the  survival  function 

#  and  its  derivatives  wrt  a  and  aa  for  the  standard 

#  gamma  distribution  with  shape  parameter  a 

#  and  vector  x.  The  input  ep  is  for  precision  control. 

len  <-  length(x) 

S  <-  1  -  pgamma(x,  a) 
vl  <-  X  <  8 
v2  <-  8  <=  X  &  X  <  aO 
v3  <-  aO  <=  X 
Sa  <-  Saa  <-  rep(0,  len) 
if(sum(vl  ==  T)  >  0)  { 
if(a  >=  aO)  { 

Sa[vl]  <-  surva.gam(a,  x[vl],  aO) 

Saa[vl]  <-  survaa.gam(a,  x[vl],  aO)} 

else  { 

rl  <-  ceiling(aO  -  a) 

Sa[vl]  <-  surva.gam(a  +  rl,  x[vl],  aO)  -  trans2(a,  x[vl],  rl) 

Saa[vl]  <-  survaa.gam(a  +  rl,  x[vl],  aO)  -  trans3(a,  x[vl],  rl)}} 
if(sum(v2  ==  T)  >  0)  { 
if(a  >=  2  *  aO)  { 

Sa[v2]  <-  surva.gam(a,  x[v2],  aO) 

Saa[v2]  <-  survaa.gam(a,  x[v2],  aO)} 

else  { 

r2  <-  ceiling(2  *  aO  -  a) 

Sa[v2]  <-  surva.gam(a  +  r2,  x[v2],  aO)  -  trans2(a,  x[v2],  r2) 

Saa[v2]  <-  survaa.gam(a  +  r2,  x[v2],  aO)  -  trans3(a,  x[v2],  r2)}} 
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if(sum(v3  ==  T)  >  0)  { 
if(a  >=  3  *  aO)  { 

Sa[v3]  <-  surva.gam(a,  x[v3],  aO) 

Saa[v3]  <-  survaa.gam(a,  x[v3],  aO)} 

else  { 

r3  <-  eeiling(3  *  aO  -  a) 

Sa[v3]  <-  surva.gam(a  +  r3,  x[v3],  aO)  -  trans2(a,  x[v3],  r3) 
Saa[v3]  <-  survaa.gam(a  +  r3,  x[v3],  aO)  -  trans3(a,  x[v3],  r3)}} 
out  <-  rbind(S,  Sa,  Saa) 
out} 

surv.gam 
funotion(a,  x,  aO) 

{#  fname  is  surv.gam 

#  ereates  first  1 6  terms  of  the  power  series  expansion  of 

#  the  tail  of  the  Ineomplete  gamma  funetion.  Useful 

#  when  a  is  large;  tmp  is  neg  of  gammastar 

ind  <-  0:15 
k  <-  length(x) 

X  <-  outer(x,  ind, 

fae  <-  matrix(gamma(a  +  1  +  ind),  k,  16,  byrow  =  T) 
tmp  <-  exp(  -  x)  *  x^a  *  apply(X/fao,  1,  sum) 

S  <-  1  -  tmp 
S} 

surva.gam 
funotion(a,  x,  aO) 

{#  fname  is  surva.gam 

#  ereates  first  1 6  terms  of  the  power  series  expansion  of 

#  the  derivative  of  the  Ineomplete  gamma  funetion.  Useful 

#  when  a  is  large;  tmp  is  neg  of  gammastar  sub  a 

ind  <-  0:15 
k  <-  length(x) 

X  <-  outer(x,  ind,  "^") 

fae  <-  matrix(psi(a  +  1  +  md)[l,  ]/gamma(a  +  1  +  ind),  k,  16,  byrow  =  T) 
tmp  <-  exp(  -  x)  *  x^a  *  apply(X  *  fae,  1,  sum) 

Sa  <-  log(x)  *  ( -  pgamma(x,  a))  +  tmp 
Sa} 

survaa.gam 
funotion(a,  x,  aO) 

{#  fname  is  survaa.gam 

#  deals  with  Saa  for  a  large.  Ineludes  the  sum  of 

#  the  first  16  terms  in  the  power  series  expansion 

#  of  X  to  the  a  times  gammastar  sub  aa 

ind  <-0;15 
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k  <-  length(x) 

X  <-  outer(x,  a  +  ind,  "^")  #  Pa  <-  exp(  -  x)  *  apply(X/matrix(gamma(a  +  1  + 

ind),  k,  16,  byrow  =  T),  1,  sum) 

brae  <-  outer(log(x),  psi(a  +  1  +  ind)[l,  ], 

den  <-  matrix(gamma(a  +  1  +  ind),  k,  16,  byrow  =  T) 

Imp  <-  matrix(psi(a  +  1  +  ind)[2,  ],  k,  16,  byrow  =  T) 

Saa  <-  exp(  -  x)  *  apply((X  *  (brac^2  -  tmp))/den,  1,  sum) 

Saa} 


transl 

function(a,  x,  r) 

{#  fname  is  transl 

#  finite  series  transition  of  standard  gamma  edf 

#  from  small  shape  parameter  to  large  shape  parameter. 

#  The  output  is  P(a-r,x)-P(a,x)=  Inc 

k  <-  length(x) 
ind  <-  1  ;r 
N  <-  length(ind) 

X  <-  outer(x,  a  -  ind,  "^") 

den  <-  matrix(gamma(a  +  1  -  ind),  k,  N,  byrow  =  T) 
S  <-  Ine  <-  exp(  -  x)  *  apply(X/den,  1 ,  sum) 

S} 


trans2 

function(a,  x,  r) 

{#  fname  is  trans2 

#  finite  series  transition  of  the  derivative  of  the  standard 

#  gamma  survivor  function  from  small  shape  parameter  to  large 

#  shape  parameter.  The  output  is  Sa(a+r,x)-Sa(a,x)=lnca 

k  <-  length(x) 
ind  <-  1  ;r 

first  <-  matrix(0,  r,  k) 
for(j  in  l:k) 

first[,  j]  <-  dgamma(x[j],  a  +  ind) 
see  <-  matrix(0,  r,  k) 
for(j  in  l:k) 

sec[,  j]  <-  dgamma(x[j],  a  +  ind)  *  psi(a  +  ind)[l,  ] 

Sa  <-  Inea  <-  log(x)  *  apply(first,  2,  sum)  -  apply(sec,  2,  sum) 
Sa} 

transS 

function(a,  x,  r) 

{#  fname  is  transS 

#  finite  series  transition  of  the  2nd  derivative  of  the  standard 

#  gamma  survivor  function  from  small  shape  parameter  to  large 

#  shape  parameter.  The  output  is  Saa(a+r,x)-Saa(a,x)=lnca 

k  <-  length(x) 
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ind  <-  1  ;r 

dens  <-  matrix(0,  k,  r) 

brae  <-  outer(log(x),  psi(a  +  ind)[l,  ], 

for(i  in  1  :k) 

dens[i,  ]  <-  dgamma(x[i],  a  +  ind) 
see  <-  matrix(psi(a  +  ind)[2,  ],  k,  r,  byrow  =  T) 

Saa  <-  Ineaa  <-  apply(dens  *  (brac^2  -  see),  1,  sum) 
Saa} 
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Appendix  E 


Listings  of  Miscellaneous  Code 


gam.est 

function(x,  ep  =  le-006) 

{#  fname  is  gam.est;  programmer:  R.  Read 

#  returns  the  mle  estimates  for  the  gamma  distribution 

#  base  upon  the  data  x.  Newton-Raphson  iteration. 

#  output  is  maximum  likelihood  for  shape,  rate,  n,  method 

#  of  moments  for  shape,  rate 

n  <-  length(x) 
xb  <-  mean(x) 

Inbx  <-  mean(log(x)) 
ss  <-  var(x) 
rate  <-  xb/ss 

lamO  <-  rate  #  initialize  with  meth  of  moments  ests. 

shape  <-  xb^2/ss 

xO  <-  Inbx  -  log(xb) 

alph  <-  shape 

repeat  { 

rO  <-  shape 
lam  <-  shape/xb 

g  <- log(shape)  +  xO  -  psi(shape)[l,  ] 
gp  <-  1/shape  -  psi(shape)[2,  ] 
shape  <-  shape  -  g/gp 
if(abs(shape  -  rO)  <  ep) 
break} 

rate  <-  shape/xb 

out  <-  o(shape,  rate,  n,  alph,  lamO) 

names(out)  <-  e("alpha-mle",  "lambda-mle",  "n",  "alpha-mm",  "lambda-mm") 
out} 


The  function  ellipse  returns  a  three-column  matrix.  Column  1  contains  the  horizontal  and 
the  other  two  columns  provide  the  two  lobes  for  the  ellipse.  The  user  must  construct  the 
plot. 

ellipse 

function(q,  m,  d,  nO  =  100) 

{#  fname  is  ellipse.  Revised  December  2002 

#  q  is  the  matrix  of  the  quadratic  form  (2x2):  (x,  y)q(x,  y)' 

#  m  is  the  centering  vector;  (m[x],  m[y]) 

#  d  is  the  (squared)  distance  value  for  the  contour 

#  nO  is  the  half-size  of  the  number  of  points  in  each  lobe 

#  of  the  ellipse. 
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a<-  q[l,l] 
b<-  q[l,2] 
cl<-  q[2,  2] 
det  <-  a  *  cl  -  b  *  b 

xO  <-  sqrt((cl  *  d)/det)  -  le-008  #  =  max(x);  latter  term  helps  insure  that  the 

fuzz  won't  stop  the  program 

X  <-  sqrt((seq(l,  nO,  1)  *  xO)/nO)  *  sqrt(xO)  #  Refines  the  partition  for  the 
extreme  values  of  x. 

X  <-  o(  -  rev(x),  0,  x) 

rad  <-  sqrt(d  *  el  -  det  *  x  *  x) 

yl  <-  ( -  b  *  x)/el  +  rad/el 

y2  <-  ( -  b  *  x)/el  -  rad/el 

X  <-  X  +  m[l] 

yl  <-  yl  +  m[2] 

y2  <-  y2  +  m[2] 

z  <-  matrix(e(x,  yl,  yl),  neol  =  3) 

z} 


The  following  funetion  produees  Figure  3.1.  Onee  exeeuted,  the  user  must  point  the 
cursor  to  a  place  in  the  lower  right  spaee  of  the  plots  and  eliek.  This  will  print  the 
legend  and  elose  out  the  funetion. 

expl.lnorm 

funotion(mu  =  0,  sig  =  1) 

{#  fname  is  expl.lnorm 

#  survival  funetions  plotted  to  eompare  Exp(l)  w/lognormal(mu,  sig) 

X  <-  o(seq(0.02,  0.55,  0.01),  seq(0.56,  3.6,  0.02)) 
y  <-  dnorm((log(x)  -  mu)/sig)/(x  *  sig) 

SI  <-  1  -  pnorm(log(x)  -  mu)/sig 
he  <-  rep(l,  length(x)) 

XX  <-  x[-(l;27)] 

Se  <-  exp(  -  xx) 
yy  <-  exp(  -  xx) 
pari  <-  par 
par(mfrow  =  e(3,  2)) 

plot(x,  y,  ylim  =  e(0,  0.75),  ylab  =  "density  funetion",  type  =  "1") 
hnes(xx,  yy,  Ity  =  4) 

title(main  =  "Compare  Density  funetions  ") 

plot(x,  exp(  -  x),  ylab  =  "survivor  funetion",  type  =  "1",  Ity  =  4) 

title(main  =  "Compare  Survivor  Funetions  ") 

hnes(x,  SI) 

plot(x,  y/Sl,  ylim  =  o(0,  1.1),  ylab  =  "hazard  funetion",  type  =  "1") 
hnes(x,  he,  Ity  =  4) 

title(main  =  "Compare  Hazard  Functions  ") 
plot(x,  X,  ylab  =  "cum  haz  function",  type  =  "1",  Ity  =  4) 
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lines(x,  -  log(Sl)) 

title(main  =  "Compare  Cum.  Hazard  Functions  ") 

SSI  <-  exp(mu  +  0.5  *  sig)  *  (1  -  pnorm(x  -  sig))  -  exp(mu  +  sig  *  x)  * 
(1  -pnorm(x)) 

LI  <-  SSl/Sl 
a  <-  max(Ll) 
bl  <-  min(x) 
b2  <-  max(x) 

plot(x,  LI,  ylim  =  c(0,  1.5  *  a),  ylab  =  "expected  resid.  life",  type  =  "1") 
lines(x,  he,  Ity  =  4) 

title(main  =  "Compare  E  {Residual  Life}  ") 

X  <-  c(0,  1,  1,  0,  0) 

y  <-  rev(x)  #  plot(box(),  xlab  =  ylab  =  axes  =  F) 

par(mar  =  c(4.9,  3,  3.8,  0.9)) 

plot(x,  y,  type  =  "1",  xlab  =  ylab  =  axes  =  F) 

lines(c(0,  0),  c(l,  0),  Ity  =  1) 

lines(c(0,  1),  c(l,  1),  Ity  =  4)  #text(0.4,  0.4,  Ftext) 

title(main  =  "Fegend") 

par  <-  pari 

leg.names  <- c("Fognormal(0,  1)  "  mean  =1.65  stdev  =  1.68", 

"Exponential}  1)  ", 

"  mean  =  1  stdev  =  1") 
legend(locator(l),  leg.names,  Ity  =  c(l,  0,  4,  0))} 


The  following  code  produces  Figure  3.2.  Once  plotted  the  user  must  click  on  the 
legend  box  in  order  to  place  information  there  and  close  out  the  function. 

weil.gam 

function} alph  =  0.75,  beta  =  2) 

{#  finame  is  weil  .gam 

#  survival  functions  plotted  to  compare  Weibull}3/4,  2)  w/gamma}.72,  rate=  .28) 

#  gamma  parameters  are  mle's  from  simulated  rweibull}200,  3/4,  2) 

X  <-  c}0.01,  seq}0.02,  15,  0.08)) 
y  <-  dweibull}x,  0.75,  2) 

Sw  <-  1  -  pweibull}x,  0.75,  2) 

Sg  <-  1  -  pgamma}x,  0.72,  0.28) 
yy  <-  dgamma}x,  0.72,  0.28) 
pari  <-  par 
par}mfrow  =  c}3,  2)) 

plot}x,  y,  ylim  =  c}0,  1.45),  ylab  =  "density  function",  type  =  "1") 
hnes}x,  yy,  Ity  =  4) 

title}main  =  "Compare  Density  functions  ") 
plot}x,  Sw,  ylab  =  "survivor  function",  type  =  "1") 
title}main  =  "Compare  Survivor  Functions  ") 
hnes}x,  Sg,  Ity  =  4) 
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plot(x,  y/Sw,  ylim  =  c(0,  1.45),  ylab  =  "hazard  function",  type  =  "1") 
lines(x,  yy/Sg,  Ity  =  4) 

title(main  =  "Compare  Hazard  Functions  ") 

plot(x,  -  log(Sw),  ylab  =  "cum  haz  function",  type  =  "1",  ylim  =  c(0,  5)) 
lines(x,  -  log(Sg),  Ity  =  4) 

title(main  =  "Compare  Cum.  Hazard  Functions  ") 

SSw  <-  (8/3)  *  gamma(4/3)  *  (1  -  pgamma((0.5  *  x)^0.75,  4/3)) 

Lw  <-  SSw/Sw 

plot(x,  Lw,  ylab  =  "expected  resid.  life",  type  =  "1") 

SSg  <-  -  X  *  Sg  +  (0.72/0.28)  *  (1  -  pgamma(x,  1.72,  0.28)) 

Lg  <-SSg/Sg 
lines(x,  Lg,  Ity  =  4) 

title(main  =  "Compare  E  {Residual  Life}  ") 

X  <-  c(0,  1,  1,  0,  0) 

y  <-  rev(x)  #  plot(box(),  xlab  =  ylab  =  axes  =  F) 

par(mar  =  c(4.9,  3,  3.8,  0.9)) 

plot(x,  y,  type  =  "1",  xlab  =  ylab  =  axes  =  F) 

lines(c(0,  0),  c(l,  0),  Ity  =  1) 

lines(e(0,  1),  c(l,  1),  Ity  =  4)  #text(0.4,  0.4,  Ltext) 

title(main  =  "Legend") 

par  <-  pari 

leg.names  <- c("Weibull(0.75,  2)  "  mean  =1.84  stdev  =  2.33", 

"Gamma(.72,  .28)  ",  "  mean  =  2.57  stdev  =  3.03") 
legend(loeator(l),  leg.names,  Ity  =  c(l,  0,  4,  0))} 


The  following  function  produces  Figure  3.3.  Onee  plotted  the  user  must  cliek  on  the 
legend  box  in  order  to  place  information  there  and  elose  out  the  funetion. 

wei2.gam 

function(alph  =  2,  beta  =  2) 

{#  finame  is  wei2.gam 

#  survival  functions  plotted  to  eompare  Weibull(2,  2)  w/gamma(3.75,  rate  =  2) 

#  gamma  parameters  are  close  to  the  mle's  from  simulated  rweibull(200,  2,  2) 

X  <-  seq(0,  5,  0.025) 
y  <-  dweibull(x,  2,  2) 

Sw  <-  1  -  pweibull(x,  2,  2) 

Sg  <-  1  -  pgamma(x,  3.75,  2) 
yy  <-  dgamma(x,  3.75,  2) 
pari  <-  par 
par(mfrow  =  e(3,  2)) 

plot(x,  y,  ylim  =  c(0,  0.48),  ylab  =  "density  function",  type  =  "1") 
lines(x,  yy,  Ity  =  4) 

title(main  =  "Compare  Density  funetions  ") 
plot(x,  Sw,  ylab  =  "survivor  function",  type  =  "1") 
title(main  =  "Compare  Survivor  Funetions  ") 
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lines(x,  Sg,  Ity  =  4) 

plot(x,  y/Sw,  ylab  =  "hazard  function",  type  =  "1") 
lines(x,  yy/Sg,  Ity  =  4) 

title(main  =  "Compare  Hazard  Funetions  ") 

plot(x,  -  log(Sw),  ylab  =  "cum  haz  function",  type  =  "1") 

lines(x,  -  log(Sg),  Ity  =  4) 

title(main  =  "Compare  Cum.  Hazard  Functions  ") 

SSw  <-  gamma(0.5)  *  (1  -  pgamma((0.5  *  x)^2,  0.5)) 

Lw  <-  SSw/Sw 
yM  <-  max(Lg) 

SSg  <-  -  X  *  Sg  +  (3.75/2)  *  (1  -  pgamma(x,  4.75,  2)) 

Lg  <-SSg/Sg 
ym  <-  min(Lw) 

plot(x,  Lw,  ylab  =  "expected  resid.  life",  type  =  "1",  ylim  =  c(ym,  yM)) 
lines(x,  Lg,  Ity  =  4) 

title(main  =  "Compare  E  {Residual  Life}  ") 

X  <-  e(0,  1,  1,  0,  0) 

y  <-  rev(x)  #  plot(box(),  xlab  =  ylab  =  axes  =  F) 

par(mar  =  e(4.9,  3,  3.8,  0.9)) 

plot(x,  y,  type  =  "1",  xlab  =  ylab  =  axes  =  F) 

lines(c(0,  0),  e(l,  0),  Ity  =  1) 

lines(e(0,  1),  c(l,  1),  Ity  =  4)  #text(0.4,  0.4,  Ltext) 

title(main  =  "Legend") 

par  <-  pari 

leg.names  <-  e("Weibull(2,  2)  ",  "  mean  =  1.77  stdev  =  2.18",  "Gamma(3.75,  2) 
", "  mean  =1.88  stdev  =  0.97") 

legend(loeator(l),  leg.names,  Ity  =  e(l,  0,  4,  0))} 

wei.eost 

funetion(rat,  alpha,  beta) 

{#  finame  is  wei.eost 

#  rat  is  the  ratio  of 

XX  <-  ceiling(log(500)) 
x2  <-  beta  *  xx^(  1/alpha) 

X  <-  seq(0,  x2,  length  =  200) 
y  <-  (x/beta)^alpha 

cost  <-  ((alpha/beta)  *  (1  +  rat  *  pgamma(y,  l)))/(gamma(l /alpha)  *  pgamma(y, 

1 /alpha)) 

m  <-  min(cost) 

M  <-  max(x) 

ind  <-  match(m,  cost)  #  print(x[ind],  m) 

#  plot(x,  cost,  type  =  "1") 

out  <-  o(ind,  M,  m) 
out} 

gam.eost 
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function(rat,  alpha,  lam) 

{#  fname  is  gam.  cost 

#  rat  is  the  ratio 

#  alpha  &  lam  are  the  parameters  of  the  gamma  distribution 

#  output  is  the  min  eost  planned  replacement  poliey 

x2  <-  alpha/lam  +  (3  *  sqrt(alpha))/lam 
X  <-  seq(0.02,  x2,  length  =  200) 

Sg  <-  1  -  pgamma(x,  alpha,  lam) 

eost  <-  (1  +  rat  *  pgamma(x,  alpha,  lam))/(x  *  Sg  +  (alpha/lam)  *  pgamma(x, 

1  +  alpha,  lam)) 
m  <-  mm(oost) 

M  <-  max(x) 

ind  <-  match(m,  eost) 

xO  <-  round(x[md],  2) 

prmt(o(ind,  xO,  m))  #  plot(x,  eost,  type  =  "1") 

out  <-  o(md,  xO,  M,  m) 

out} 

NormCT 

funotion(x,  alph,  graph  =  T) 

{#  fname  is  NormCT 

#  this  is  a  plotting  funetion  that  produees  1-alph  level 

#  eonfidenee  trapezoids  for  (mu,  sig)  for  normal  data  x. 

#  If  graph  =  F  then  the  output  is  a  row  matrix  eontaining 

#  the  vertiees  of  the  trapezoid. 

out  <-  NULL 

n  <-  length(x) 

alphl  <-  1  -  sqrt(l  -  alph) 

zO  <-  -  qnorm( alph  1/2) 

kl  <-  qohisq(alphl/2,  n  -  1) 

k2  <-  qchisq(l  -  alphl/2,  n  -  1) 

SS  <-  (n  -  1)  *  var(x) 

xb  <-  mean(x) 

s  <-  stdev(x) 

sigl  <-  sqrt(SS/k2) 

sigu  <-  sqrt(SS/kl) 

xlmin  <-  xb  -  (zO  *  sigl)/sqrt(n) 

xlmax  <-  xb  +  (zO  *  sigl)/sqrt(n) 

xumin  <-  xb  -  (zO  *  sigu)/sqrt(n) 

xumax  <-  xb  +  (zO  *  sigu)/sqrt(n) 

XX  <-  o(xlmin,  xlmax,  xumax,  xumin,  xlmin) 
yy  <-  o(sigl,  sigl,  sigu,  sigu,  sigl) 
if(graph  ==  T)  { 

plot(xx,  yy,  type  =  "1",  xlab  =  "mu",  ylab  =  "sigma") 
title(mam  =  "Confidenee  Trapezoid  for  Normal  Distribution") 
points(xb,  s)} 
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else  out  <-  rbind(xx,  yy) 
return(out)} 


area.eomp 

funetion(x,  w,  uO  =  0,  yO  =  0,  rnd  =  4) 

{#  fname  is  area.eomp 

#  Computes  the  signed  net  areas  separating  the  empirieal 

#  edf  s  of  the  ordered  sets  x  and  w.  These  edf  s  are  polygonal 

#  curves  which  are  connected  with  straight  line  segments.  The 

#  two  data  sets  are  of  the  same  length.  It  seems  necessary  to 

#  do  some  rounding  because  of  "fuzz"  problems;  hence  the  input 

#  rounding  quantity,  "rnd". 

nO  <-  length(x) 

X  <-  round(x,  md) 

w  <-  round(w,  md) 

out  <-  matrix(0,  ceilmg((l  +  n0)/2),  5) 

jj  <-  1 

repeat  { 

out[jj,  ]  <-  seg.comp(x,  w,  uO,  yO,  nO,  jj) 
uO  <-  round(out[jj,  4],  md) 
tmp  <-  round(out[jj,  5],  rnd) 
yO  <-  tmp  -  floor(tmp) 

X  <-  x[x  >  uO] 
w  <-  w[w  >  uO] 
nl  <-  length(x) 
n2  <-  length(w) 
if(nl  ==  0  I  n2  ==  0) 
break 

if(all(x  ==  w)) 
break 

jj  <-jj  +  1} 

out<- out[l;jj,  ] 
out} 


seg.comp 

function(x,  w,  uO,  yO,  nO,  jj) 

{#  fname  is  seg.comp 

#  Computes  the  areas  under  the  polygonal  curves,  between 

#  two  knots,  and  returns  their  difference  (signed).  A  flag 

#  is  set  =  0  if  the  x  edf  is  above  the  w  edf,  and  set  =  1 

#  otherwise.  The  x  and  w  vectors  are  mono  increasing;  n  is 

#  the  number  of  points  in  the  full  sets.  The  initial  points 

#  (uO,  yO)  mark  the  beginning  of  the  segment;  the  crossover 

#  point  (ul,  yl)  is  the  segment  end  and  is  computed  internally; 
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#  when  segment  is  open  ended,  it  is  (x[n],  n).  The  marker  fl 

#  is  set  to  one  when  the  segment  is  elose  on  the  right. 

#  f2  =  1  means  the  right  knot  is  not  a  data  point. 

n  <-  f  <-  yl  <-  length(x)  #  initialization 
fl  <-  flag  <-  rect  <-  f2  <-  0 
ul  <-  max(x[l],  w[l]) 

ss  <-  sort(c(x,  w))  #  first  look  for  superpositions  and  remove 

repeat  { 

if(x[l]!=w[l]) 
break 
uO  <-  x[l] 
x<-  x[-l] 
w  <-  w[-l] 
ss  <-  ss[  -  (1:2)] 
n  <-  f  <-  n  -  1 
if(n  ==  1  I  n  ==  0)  { 
f<-  1 
yl  <-  1 
break} 

else  { 

if(x[l]==w[l]){ 

ss  <-  ss[  -  (1:2)] 
uO  <-  x[l] 
x<-  x[-l] 
w  <-  w[-l] 

n  <-  f  <-  yl  <-  length(x) 
if(n  ==  0) 
break  } } } 

#  Set  the  flag  &  initiate 

if(n  >  0) 

j  <-  1  :n 

if(length(ss)  >  0)  { 

ind  <-j[x|j]  >=w|j]] 
if(w[l]==ss[l]){ 
flag  <-  1 

ind<-j[w|j]  >=x|j]]  } 
if(is.na(ind[l])  ==  F)  { 
fl  <-  1 
f<-  ind[l]}} 

#  initialize  the  end  eorreetions  and  the  eenter  eomputations 

areal  <-  area2  <-  adj  1  <-  adj2  <-  0 
ul  <-  max(x[f],  w[f]) 
yl  <-f 

if(f  >  1  &  fl  ==  1)  { 

Pl<-  c(x[(f-l):f]) 

P2  <-  e(w[(f-  l):f]) 
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tout  <-  sol.pt(Pl,  P2) 
ul  <-  tout[l] 
yl  <-  tout[2]  +  f 
if(x[f]  !=w[f])  { 

£2<-  1 
f<-f-l}} 

#  prepare  for  the  open  ended  ease 

if(n  >  0)  {reet  <-  f  *  abs(x[n]  -  w[n]) 

j  <-  1  :n 

areal  <-  ((x[l]  -  uO)  *  (1  +  y0))/2 
area2  <-  ((w[l]  -  uO)  *  (1  +  y0))/2 

if(f  >  1)  {areal  <-  areal  +  ((fl  *  f2  *  (ul  -  x[f  -  1])  *  (yl  +  f  -  l))/2) 

area2  <-  area2  +  ((fl  *  f2  *  (ul  -  w[f  -  1])  *  (yl  +  f  -  l))/2) 

ff  <-  f  #  prep  adiustment  for  existenee  of  interior  part  of  segment 
if(x[f]  !=w[f]) 

ff<-ff-  1 

adjl  <-  0.5  *  (x[ff]  *  (2  *  ff  +  1)  -  x[l])  -  sum(x[l:ff]) 
adj2  <-  0.5  *  (w[ff]  *  (2  *  ff  +  1)  -  w[l])  -  sum(w[l:ff])  } 
areal  <-  areal  +  (1  -  fl)  *  (1  -  flag)  *  reet 
area2  <-  area2  +  (1  -  fl)  *  flag  *  reet  } 
net  <-  (areal  +  adj  1  -  area2  -  adj2)/n0 
out  <-  o(net,  flag,  f,  ul,  yl) 

names(out)  <-  o("net",  "flag",  "seg. length",  "hend",  "vend") 
out} 


sol.pt 

funotion(Pl,  P2) 

{#  finame  is  sol.pt 

#  finds  the  erossover  solution  point 

#  for  two  edf  s  that  have  the  same  number 

#  of  pts  in  the  horiz  &  equi  spaeed  in  the  vertieal. 

#  If  both  delx  and  delq  are  zero,  then  y  is  set  to 

#  y  =  -1 .  Otherwise  0  <  =  y  <  1 

xl  <-Pl[l] 

x2  <-Pl[2] 

wl  <-  P2[l] 

w2  <-  P2[2] 

delx  <-  x2  -  xl 

delw  <-  w2  -  wl 

if(delx  ==  0  &  delw  ==  0)  { 

X  <-  xl 
y<--l  } 

if(delx  ==  0  &  delw  >  0)  { 

X  <-  xl 
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y  <-  (x  -  wiydelw  } 
if(delx  >  0  &  delw  ==  0)  { 

X  <-  wl 

y  <-  (x  -  xiydelx  } 
if(delx  >  0  &  delw  >  0)  { 

X  <-  (xl  *  delw  -  wl  *  delxy(delw  -  delx) 
y  <-  (x  -  wiydelw 
if(xl  ==  wl) 

y<-0  } 

out  <-  e(x,  y) 
out} 

sol.pt 

funotion(Pl,  P2) 

{#  fname  is  sol.pt 

#  finds  the  erossover  solution  point 

#  for  two  edf  s  that  have  the  same  number 

#  of  pts  in  the  horiz  &  equi  spaeed  in  the  vertieal. 

#  If  both  delx  and  delq  are  zero,  then  y  is  set  to 

#  y  =  -1 .  Otherwise  0  <  =  y  <  1 

xl  <-Pl[l] 

x2  <-Pl[2] 

wl  <-  P2[l] 

w2  <-  P2[2] 

delx  <-  x2  -  xl 

delw  <-  w2  -  wl 

if(delx  ==  0  &  delw  ==  0)  { 

X  <-  xl 
y<--l  } 

if(delx  ==  0  &  delw  >  0)  { 

X  <-  xl 

y  <-  (x  -  wl)/delw} 
if(delx  >  0  &  delw  ==  0)  { 

X  <-  wl 

y  <-  (x  -  xl)/delx} 
if(delx  >  0  &  delw  >  0)  { 

X  <-  (xl  *  delw  -  wl  *  delx)/(delw  -  delx) 
y  <-  (x  -  wl)/delw 
if(xl  ==  wl) 
y<-0} 

out  <-  o(x,  y) 
out} 

Newt.wei 

funotion(x,  t,  param,  ep  =  0.0001) 

{#  fname  is  Newt.wei 
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#  Newton-Raphson  method,  then  biseetion  seareh  used  to  find  mle 

#  for  the  Weibull  distribution  using  the  extreme  value  distribution  teehnique. 

#  Data  input  is  (x,  t);  param  initialization  is  (u,  b),  b  >  0 

#  output  is  (u,  b,  alph,  beta),  r  is  eardinality  of  uneensored  set 

r  <-  length(x) 

V  <-  log(x) 
vb  <-  mean(v) 

It  <-  log(t) 
b  <-  param[2] 
u  <-  param[l] 
vv  <-  e(v,  It) 
b2  <-bl  <-b 
hi  <-  h2  <-  0 

j<-l 

flag  <-  1 
repeat  { 

if(flag  ==  1)  { 

D  <-  sum(exp(vv/b)) 

N  <-  sum(vv^2  *  exp(vv/b)) 

A  <-  sum(vv  *  exp(vv/b)) 
h  <-  A/D  -  b  -  vb 

hp  <-  -  N/(D  *  b^2)  +  A^2/(D  *  b)^2  -  1 
b<-bl  -h/hp 
if(b  <=  0) 

b<-0.01 

if(max(abs(h),  abs(b  -  bl))  <  ep) 
break 

if(j/2  !=  round(j/2))  { 
bl  <-b 
hi  <-h} 

if(j/2  ==  round(j/2))  { 
b2  <-b 
h2<-h} 

j  <-  j  +  1 

if(sign(h2  *  hi)  <  0) 
flag  <-  0} 

ifG==25) 

break 
if(flag  ==  0)  { 

b  <-  (bl  +  b2)/2 
D  <-  sum(exp(vv/b)) 

A  <-  sum(vv  *  exp(vv/b)) 
h  <-  A/D  -  b  -  vb 
if(sign(h  *  hi  <  0))  { 
h2  <-h 
b2  <-b} 
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if(sign(h  *  h2  <  0))  { 
hi  <-h 
bl  <-b} 

if(max(abs(h),  abs(b  -  bl))  <  ep) 
break 

j  <-j  +  1 

ifO  ==  50) 

break}  } 

u  <-  b  *  log(D/r) 
alph  <-  l^ 
beta  <-  exp(u) 

out  <-  e(u,  b,  alph,  beta,  flag) 

names(out)  <-  e("u",  "b",  "alph",  "beta",  "flag") 

return(out)} 


exp.lnorm 

funetion(mu  =  0,  sig  =  1) 

{#  fname  is  exp.lnorm 

#  survival  funetions  plotted  to  eompare  Exp(l)  w/lognormal(mu,  sig) 
X  <-  o(seq(0.02,  0.55,  0.01),  seq(0.56,  3.6,  0.02)) 
y  <-  dnorm((log(x)  -  mu)/sig)/(x  *  sig) 

SI  <-  1  -  pnorm(log(x)  -  mu)/sig 
he  <-  rep(l,  length(x)) 

XX  <-  x[-(l;27)] 

Se  <-  exp(  -  xx) 
yy  <-  exp(  -  xx) 
split.sereen(e(3,  2)) 
oldpar  <-  par() 
on.exit(par(oldpar)) 
sereen(l) 

par(eex  =  0.8,  mar  =  e(5,  6,  4,  2)) 

plot(x,  y,  ylim  =  e(0,  0.75),  ylab  =  "density  funetion",  type  =  "1") 
lines(xx,  yy,  Ity  =  4) 

title(main  =  "Compare  Density  funetions  ") 
soreen(2) 

par(oex  =  0.8,  mar  =  o(5,  6,  4,  2)) 

plot(x,  exp(  -  x),  ylab  =  "survivor  funetion",  type  =  "1",  Ity  =  4) 
title(main  =  "Compare  Survivor  Functions  ") 
lines(x,  SI) 
screen(3) 

par(cex  =  0.8,  mar  =  c(5,  6,  4,  2)) 

plot(x,  y/Sl,  ylim  =  c(0,  1.1),  ylab  =  "hazard  function",  type  =  "1") 
lines(x,  he,  Ity  =  4) 

title(main  =  "Compare  Hazard  Functions  ") 
screen(4) 
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par(cex  =  0.8,  mar  =  c(5,  6,  4,  2)) 

plot(x,  X,  ylab  =  "cum  haz  function",  type  =  "1",  Ity  =  4) 

lmes(x,  -  log(Sl)) 

title(mam  =  "Compare  Cum.  Hazard  Functions  ") 
screen(5) 

par(cex  =  0.8,  mar  =  c(5,  6,  4,  2)) 

SSI  <-  exp(mu  +  0.5  *  sig)  *  (1  -  pnorm(x  -  sig))  -  exp(mu  +  sig  *  x)  * 

(1  -pnorm(x)) 

LI  <-  SSl/Sl 
a  <-  max(Ll) 
bl  <-  min(x) 
b2  <-  max(x) 

plot(x,  LI,  ylim  =  c(0,  1.5  *  a),  ylab  =  "expected  resid.  life",  type  =  "1") 
lines(x,  he,  Ity  =  4) 

title(main  =  "Compare  E  {Residual  Life}  ") 
screen(6) 

par(cex  =  0.8,  mar  =  c(5,  6,  4,  2)) 

X  <-  c(0,  1,  1,  0,  0) 

y  <-  rev(x)  #  plot(box(),  xlab  =  ylab  =  axes  =  F) 
plot(x,  y,  type  =  "n",  xlab  =  ylab  =  axes  =  F)  ##  lines(c(0,  0),  c(l,  0), 
Ity  =  1) 

###  lines(c(0,  1),  c(l,  1),  Ity  =  4)  #text(0.4,  0.4,  Ftext) 

leg.names  <-  c("Fognormal(0,  1)",  "mean  =  1.65",  "sd  =  1.68",  "Exponential}!)", 
"mean=  1",  "sd=  1",  "",  "") 
leg.xy  <-  locator}  1) 

legend}leg.xy$x,  leg.xySy,  leg.names,  Ity  =  c}l,  0,  0,  4,  0,  0,  0,  0),  cex  =1) 
close. screen}all  =  T)} 

ext.newt 

function} input,  dat,  ep  =  10^-4) 

{#  finame  is  ext.newt 

#  Newton's  method  applied  to  extreme  value  distribution 
u  <-  uO  <-  input}  1] 
b  <-  bO  <-  input[2] 

X  <-  dat[dat[,  2]  ==  1,  1] 
a  <-  dat[,  2] 

w  <-  a  *  dat[,  1]  +  }1  -  a)  *  dat[,  3] 
r  <-  sum} a) 
gO  <-0 

repeat  {g  <-  sum}w  *  exp}w/b))/sum}exp}w/b))  -  b  -  mean}x) 
gl  <-  sum}w^2  *  exp}w/b))/sum}exp}w/b)) 
g2  <-  sum}w  *  exp}w/b))/sum}exp}w/b)) 
gp<-  -  }gl  +  g2^2  -  l)/b^2 
b  <-  bO  -  g/gp 
if}b  <  0) 

b  <-  0.2 
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if(max(abs(b  -  bO),  abs(g  -  gO))  <  ep) 
break 

gO  <-g 
bO  <-b} 

u  <-  b  *  log(sum(exp(w/b))/r) 
retum(c(u,  b))} 

cost.comp 

function(rat) 

{#  fname  is  cost.comp 

#  Compares  the  cost  of  planned  replacement  curves  of  our  two  competing  IFR  models 
Weibull(2,  2)  &  gamma(3.75,  rate  =  2) 

#  Three  sets  of  curves  are  generated,  one  for  each  of  the  n  ratios  in  the  input  rat 

n  <-  length(rat) 
aw  <-  2 
bw  <-  2 
XX  <-  log(500) 

x2  <-  0.25  *  bw  *  XX  *  (1/aw) 
ag  <-  3.75 
lam  <-  2 

xl  <-  ag/lam  +  (1  *  sqrt(ag))/lam 

x3  <-  max(xl,  x2)  #  range  of  variate  values 

X  <-  seq(0,  x3,  length  =  200) 

y  <-  (x/bw)^aw  #y  <-  dweibull(x,  2,  2) 

Sw  <-  1  -  pweibull(x,  2,  2) 

Sg  <-  1  -  pgamma(x,  3.75,  2)  #yy  <-  dgamma(x,  3.75,  2) 
pari  <-  par 
par(mfrow  =  c(n,  1)) 

numw  <-  (aw/bw)  *  (1  +  outer(rat,  pgamma(y,  1),  "*")) 

denw  <-  matrix(gamma(l/aw)  *  pgamma(y,  1/aw),  byrow  =  T,  mow  =  n,  ncol  = 
200) 

costw  <-  numw/denw 

numg  <-  1  +  outer(rat,  pgamma(x,  ag,  lam),  "*") 

deng  <-  matrix(x  *  Sg  +  (ag/lam)  *  pgamma(x,  1  +  ag,  lam),  byrow  =  T,  mow  = 
n,  ncol  =  200) 

costg  <-  numg/deng 
for(j  in  l:n)  { 

plot(x,  costw[j,  ],  ylab  =  c("rel.  cost"),  xlab  =  c("time  in  service"),  type  = 

"1") 

lines(x,  costgQ,  ],  Ity  =  4) 
title(main  =  paste("ratio  =  ",  rat[j]))  } 
mcostw  <-  apply(costw,  1 ,  min) 
mcostg  <-  apply(costg,  1,  min) 
out  <-  rbind(mcostw,  mcostg)#x  <-  c(0,  1,  1,  0,  0) 

#  y  <-  rev(x)  #  plot(box(),  xlab  =  ylab  =  axes  =  F) 

#  par(mar  =  c(4.9,  3,  3.8,  0.9)) 
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#  plot(x,  y,  type  =  "1",  xlab  =  ylab  =  axes  =  F) 

#  lines(e(0,  0),  e(l,  0),  Ity  =  1) 

#  lines(e(0,  1),  e(l,  1),  Ity  =  4)  #text(0.4,  0.4,  Ltext) 

#  title(main  =  "Legend") 

par  <- pari  #leg.names  <- e("Weibull(2,  2) "  mean  =1.77  stdev  =  2.18", 

#  "Gamma(3.75,  2)  ",  "  mean  =1.88  stdev  =  0.97") 

#  legend(loeator(l),  leg. names,  Ity  =  e(l,  0,  4,  0)) 

out} 

gam.eost 

funetion(rat,  alpha,  lam) 

{#  fname  is  gam.eost 

#  rat  is  the  ratio 

#  alpha  &  lam  are  the  parameters  of  the  gamma  distribution 

#  output  is  the  min  eost  planned  replaeement  poliey 

x2  <-  alpha/lam  +  (3  *  sqrt(alpha))/lam 
X  <-  seq(0.02,  x2,  length  =  200) 

Sg  <-  1  -  pgamma(x,  alpha,  lam) 

eost  <-  (1  +  rat  *  pgamma(x,  alpha,  lam))/(x  *  Sg  +  (alpha/lam)  *  pgamma(x, 

1  +  alpha,  lam)) 
m  <-  min(oost) 

M  <-  max(x) 

ind  <-  matoh(m,  eost) 

xO  <-  round(x[ind],  2) 

print(o(ind,  xO,  m))  #plot(x,  eost,  type  =  "1") 

out  <-  o(ind,  xO,  M,  m) 

out} 

ks.dist 

funotion(dat,  mod) 

{#  fname  is  ks.dist 

#  the  distanee  between  the  data  edf  and  the  model  edf 

#  using  the  Kolmorgoroe-Smirnov  distanee  funetion.  The 

#  values  of  mod  must  be  the  edf  of  the  model  evaluated 

#  at  the  values  of  dat. 

n  <-  length(dat) 
dat  <-  sort(dat) 
mod  <-  sort(mod) 
p  <-  (l;n)/n 
pi  <-  o(0,  p[  -  n]) 

ks  <-  max(abs(o(mod  -  p,  pi  -  mod))) 
return(ks)} 

Gsq.norm 

funotion(x,  y,  xb,  s,  n) 

{#  fname  is  Gsq.norm 
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#  This  program  computes  the  vertieal  component,  z,  of  the 

#  G-squared  statistic  over  the  grid  work  that  results 

#  from  the  x,  y  values,  using  the  normal  likelihood 

#  function.  Thus,  one  can  construct  contour  plots  of 

#  (x,  y,  z).  xb  and  s  are  the  mle's  of  mu  and  sigma  from 

#  a  random  sample  of  size  n.  x  values  go  with  mu; 

#  y  with  sigma 

SS  <-  n  *  s^2 

tmpl  <-  2  *  n  *  log(y/s)  +  SS/y^2  -  n 
tmp2  <-  n  *  outer((xb  -  x)^2,  y^2,  "/") 
zz  <-  tmp  1  +  t(tmp2) 
z  <-  t(zz) 
retum(z)} 


Gsq.gam 

function(dat,  x,  y,  ah,  bh,  n) 

{#  finame  is  Gsq.gam 

#  This  program  eomputes  the  vertieal  component,  z,  of  the 

#  G-squared  statistie  over  the  grid  work  that  results 

#  from  the  x,  y  values,  using  the  gamma( alpha,  beta)likelihood 

#  funetion.  Thus,  one  ean  eonstruct  eontour  plots  of 

#  (x,y,z).  ah  and  bh  are  the  mle's  of  alpha  and  beta  from 

#  a  random  sample,  dat,  of  size  n.  x  values  go  with  alpha; 

#  y  with  beta 

xb  <-  mean(dat) 
blx  <-  mean(log(dat)) 

za  <-  -  Igamma(ah)  -  ah  *  log(bh)  +  ah  *  blx  -  ah  #  scalar  term 
zl  <-  Igamma(x)  -  x  *  blx  #  vector  terms 
z21  <-  outer(x,  log(y),  "*") 

z22  <-  xb  *  matrix(l/y,  length(x),  length(y),  byrow  =  T) 
zz  <-  za  +  zl  +  z21  +  z22 
out  <-  2  *  n  *  zz 
return(out)} 

Gsq.wei 

funotion(dat,  x,  y,  ah,  bh,  n) 

{#  finame  is  Gsq.wei 

#  This  program  computes  the  vertical  component,  z,  of  the 

#  G-squared  statistie  over  the  grid  work  that  results 

#  from  the  x,  y  values,  using  the  Weibull  likelihood 

#  function.  Thus,  one  can  construct  contour  plots  of 

#  (x,  y,  z).  ah  and  bh  are  the  mle's  of  shape  and  seale 

#  a  random  sample,  dat,  of  size  n.  x  values  go  with  ah; 

#  y  with  bh 

blx  <-  mean(log(dat)) 
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zO  <-  log(ah)  -  ah  *  log(bh)  +  ah  *  blx  -  1 

zl  <-  log(x)  +  X  *  blx  #  vectors 

V  <-  outer(x,  log(y),  "*") 

xba  <-  apply(outer(dat,  x,  2,  sum)/n 

zz  <-  zO  -  zl  +  V  +  xba  *  exp(  -  v) 

z  <-  2  *  n  *  zz 

retum(z)} 


#  scalars 
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